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Background 

This research grant topic addressed the capability of using blade-mounted 
instrumentation, notably accelerometers, to reconstruct blade modal motion and forced 
response from both wind tunnel and flight test data. Two documents produced as a 
consequence of this investigation are attached as an Appendix: the first, Reconstruc 10 
of Helicopter Rotor Blade Modal Properties and Examination of Blade-Mounted Sensor 
Data”, comprised part of a visiting student’s thesis for the degree of Diplom-Ingemeur 
from the University of Stuttgart; the second, “Instrumented Blade Experiments Using a 
Light Autogiro”, is a conference paper from the 16"’ European Rotorcraft Forum, held m 
Glasgow, Scotland in September of 1990, and is a condensed version of a Master s Thesis 
submitted to the Department of Mechanical and Aerospace Engineering at Princeton 

University. 

The first document explores, after some tutorial material, the possibility of 
identifying the forced modal response of the rotorblade using only accelerometer 
measurements. Through a collection of various filtering concepts, it is shown that some 
success may be achieved in this regard, although the method involves considerable 
“tuning” on the part of the engineer and hence is not yet suitable for “production wor c. 
A finite element analysis is developed that provides some technical basis for design of e 

various filtering methods. 

The second document discusses the developmental work associated with 
instrumenting a light autogiro's rotor blades with accelerometers for measurement of in- 
plane and out-of-plane response. Issues related to data formatting, power conditioning, 
and impact testing of the data system and sensors are discussed. The ultimate goal of the 
project - the demonstration testing of the system on the autogiro in towed operation - as 
yet to transpire, as the runway that was to serve as the test track was destroyed by the 
University to support further development of the Forrestal Campus. 

While these two works did not directly generate a solution to this complex 
identification and measurement problem, they have aided the further development of this 
concept. On-going Navy-sponsored work related to rotor instrumentation and data 
processing has been able to utilize portions of each of these works to identify methods for 
extracting both rotor motion and applied loads from similarly instrumented rotorcraft. 
Results of this work have been recently reported at the American Helicopter Societ> s 

Annual Forum. 
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Introduction 


Excessive vibration remains one one of the most difficult problems that 
faces the helicopter industry today, affecting all production helicopters at some 
phase of their development. Vibrations in rotating structures may arise from 
external periodic dynamic airloads whose frequencies are are close to the natural 
frequencies of the rotating system itself. The goal for the structures engineer 
would thus be to design a structure as free from resonance effects as possible. In 
the case of a helicopter rotor blade these dynamic loads are a consequence of 
asymmetric airload distribution on the rotor blade in forward flight, leading to a 
rich collection of higher harmonic airloads that force rotor and airframe 
response. Accurate prediction of the dynamic characteristics of a helicopter rotor 
blade will provide the opportunity to affect in a positive manner noise intensity, 
vibration level, durability, reliability and operating costs by reducing 
objectionable frequencies or moving them to a different frequency range and 
thus providing us with a lower vibration rotor. In fact, the dynamic 
characteristics tend to define the operating limits of a rotorcraft. As computing 
power has increased greatly over the last decade, researchers and engineers have 
turned to analyzing the vibrational characteristics of aerospace structures at the 
design and development stage of the production of an aircraft. Modern rotor 
blade construction methods lead to products with low mass and low inherent 
damping so careful design and analysis is required to avoid resonance and an 
undesirable dynamic performance. In addition, accurate modal analysis is 
necessary for several current approaches in elastic system identification and 
active control. 


1. Analysis of Structural Vibration 

The analysis procedure falls into the two categories: experimental 
structural analysis which involves the testing and measurement of an actually 
built structure, and analytical structural analysis, which consists of 
mathematically modeling a structure on a computer and predicting its response 
to excitations as well as inferring the inherent vibrational characteristics like 
mode shape and natural frequency. Structural analysis based on a computational 
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method allows the freedom to either experiment with a new design or change an 
existing design and then to make a statement about the magnitude of the 
changes in dynamic properties resulting from some perturbation of the original 
design. This relieves the designer from making costly physical changes to the 
system itself until he is satisfied with the characteristics presented to him in the 
simulation. One of the major problems to date is that none of the analytic 
models developed for the simulation of rotor blade behavior are capable of an 
accurate response prediction over the entire flight envelope of a helicopter. An 
exact model would have to include detailed unsteady aerodynamic predictions 
for the airloads as well as nonlinear geometric coupling effects due to the blade 
rotation. 

In response to this discrepancy between predicted and measured 
vibrational behavior, the US Army in conjunction with NASA, has undertaken 
a program of flight testing a highly instrumented UH-60A "Blackhawk" 
helicopter. The purpose of these experiments is to provide the industry and the 
university research community with a wealth of flight test data obtained under 
specific flight conditions. This data would serve to be correlated and compared 
with various parts of currently ongoing aeroelastic and vibrational helicopter 
rotor and airframe simulations. It is hoped that finding and evaluating the 
discrepancies between the database and simulation results will improve the 
understanding into where the engineering approximations and assumptions in 
the simulation programs are deficient. The database was further extended by 
conducting a static, non-rotating shake test of the same UH-60A rotor blade. In 
this shake test it was possible to determine the non-rotating dynamic properties 
of the rotor blade by identifying its modes of vibration and natural frequencies. 
As any vibrational state can be said to consist of a specific combination of its 
vibrational modes, each having a particular natural frequency, we find the area 
of coincidence of the analytical and experimental part of structure dynamics here. 
It is always important to have one set of measurements to quantify the accuracy 
of the mathematical model. The simulation model can then be adjusted or 
altered to better conform to the real model. Having once learned that the 
mathematical model does indeed correspond with the experimental data, the 
engineer can change the structural and material properties of his mathematical 
model with a high degree of confidence that his simulation calculations will 
reflect the real performance of the structure he is designing. 
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1.1 Introduction to Modal Theory 


The dynamic properties of an elastic structure are most commonly 
modeled using either a set of linear differential equations in the time domain, or 
by a set of algebraic equations after a Fourier transform of these differential 
equations of motion. The latter is commonly used during the actual testing 
phase of a design. First we will examine the former. 


1.1.1 Modal Theory 

In the analysis of dynamic systems having one or more degrees of freedom and 
their response to any sort of excitation, the theory of modal analysis can be 
applied to the structure. 

The modal theory simply states that any shape of the oscillating structure 
can be expressed in terms of natural modes of the structure and the degree of the 
participation of these modes. The natural modes are the physical shape that the 
structure takes when excited at a natural frequency. An easy way of thinking of 
these natural modes is to consider the following two degree of freedom system 
(Fig 1.1) composed of two masses that are connected by two springs and two 
damping devices. The number of natural modes that one can expect to find 
always corresponds to the number of degrees of freedom. It follows that we can 
expect to find 2 normal modes and their natural frequencies. 



Fig 1 .1 A dynamic System 
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As this structure is excited with a force, F(t), of variable frequency, we can 
envision two cases in which the masses vibrate in a stationary way. In one case, 
the masses will vibrate in phase with each other, though with different 
amplitudes. In the other case, they will vibrate 180 degrees out of phase with each 
other. Both these cases are associated with a natural frequency of the system. The 
system is a lumped parameter model, consisting of two discrete mass points. After 
deriving all the necessary equations for the lumped parameter model, we will 
expand the theory to encompass structures with distributed mass, i.e. beams. 

When we set up the equations of motion for a linear system, we obtain the 
familiar equation in matrix-form : 

Mx+Dx+Kx = F(t) (l.l.l-l) 

This equation represents the force balance between the inertial (M x) , dissipative 

(D x) and restoring forces(K x) and all outside forces(F(t)). For any model of a 
structure or system, with any number of degrees of freedom (DOF) , as long as we 
use only masses, springs and proportional dampers, we will obtain the above 
form of the equations of motion. All linear properties of a structure can be 
defined by the mass, damping and stiffness matrices. We now have a system of 
second order coupled differential equations that are linear and time-invariant. 
M , D and K are matrices of the dimension ( n x n), n being the number of DOF. 
The equation is coupled because the D and K matrices are not diagonal, only 
symmetric. To uncouple these equations we must find a basis for this system, in 
which all the matrices become diagonal, thus reducing the n coupled equations 
to a system of n uncoupled equations. This can be thought of as reducing the 
structure to a system of masses with one DOF each. The process is called modal 
decomposition. Modal decomposition in the case of light damping or if the 
damping matrix is proportional to either the mass or stiffness matrix, or 
proportional to both leads to real modes, while any damping that does not fall 
into the beforementioned categories will lead to complex modes. 

To execute such a modal decomposition we will have to find the eigenvalues, k, 
(or natural frequencies co n =^ n 0 - 5 ) and the eigenvectors, u, of our structure. In the 
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interest of simplicity we will disregard damping in the following derivations and 
solve the homogeneous form of (l.l.l-l). Assuming for x(t) 

x(t) = u e Xt (1.1. 1-2) 

and thus 

x = l 2 ue^ = ), 2 x 

Substituting into (l.l.l-l) and solving for X leads to 

a 2 M + K)u = 0 (1.1. 1-3) 

This gives us the characteristic equation (CE) for the structure. A non-trivial 
solution exists only when 

det (X 2 M + 3C)= 0 (1.1.1-4) 

This CE is a polynomial whose roots are the eigenvalues of the system. 
Depending on the kind of polynomial that the CE represents, we can imagine the 
roots to be of different character. If we assume zero damping in our example case 
we would find 4 complex roots or a pair of conjugate complex roots. The roots 
are purely imaginary, having no real part. Successively inserting these roots back 
into equation (1.1. 1-3), we receive the matrix equations for the eigenvectors u. 
These eigenvectors are only defined in terms of their direction, not their 
magnitude, since the matrix equations are linearly dependent. (Traditionally the 
eigenvectors are scaled so that some DOF that is in some way predominant is 
unity. In the case of vibrating beams we will always set that DOF to unity that 
displays the maximum deflection). The vectors thus calculated are the 
normalized eigenvectors, u n , associated with the natural frequency, w n , or the n- 
th eigenvalue, X n . In our example, we would get two of these eigenvectors, and 
could assemble them into the modal system matrix U = [ u.i , u 2 J. 1 


^Let us consider the system matrix that was found for our example: First of all we must remember 
that the eigenvectors were arbitrarily scaled, so that they represent only the amplitude 
relations of the oscillations of the lumped parameter masses. Recognizing this, we find that for 
m l =m 2' di=d}=0 and kj=k 2 that our eigenvectors are: 
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Another very helpful way of understanding the significance of eigenvectors can 
be found in Ref l.-l. The eigenvalue /eigenvector form gives an interesting 
insight into the nature of the solution: 

u , = K-l M & a-1-1-4) 

The right-hand side of this equation subjects the displacement vector at 
each lumped mass (u n ) to a matrix multiplication which should both extend and 

rotate this vector. The same vector, on the left side of the equation, is simply 
subjected to a multiplication with a scalar constant, resulting in only a expansion 
of the vector but no change in its orientation. So, the vectors we are seeking, are 
vectors that are not rotated by the matrix multiplication on the right-hand side. 
Such solution vectors are the eigenvectors. 


Using these eigenvectors, we are in the position to modaly decompose the 
equations of motion. We have said that every form that the oscillating structure 
takes can be expressed in terms of the eigenvectors : 


x(t) = U] qi(t) + M2 q2(t) = [ uj U2 ] 



9 


( 1 . 1 . 1 - 5 ) 


note: q(t) := q 

q | and q 2 give the amplitudes of the specific mode at a point in time and thus 
are responsible for the time-dependency of the equation of motion. They are 
unknown until some initial conditions are introduced to the problem. 


m= 



m= 


er) 


We see, that for the first resonant, or natural frequency, the masses are oscillating in phase 
(both components of the eigenvectors have the same sign). For the second resonant frequency the 
masses are oscillating 180 ° out of phase. 
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It can be shown that a front and back multiplication with the modal matrix, U, 
the mass and stiffness matrices are diagonalized (this uses the orthogonality 
condition, explained later). 


Mdiag= uTM U 
Kdiag = U T K U 


(l.l.l-6a) 

(l.l.l-6b) 


So, taking our original equation (1.1. 1-1) and inserting (1.1. 1-2) and then 
multiplying both sides by U we get: 


U*M U | + U t KXJq = U* F(t) 


(1.1.1 -7) 


or 


Mdiag I + Kdiag 3 = * 


( 1 . 1 . 1 - 8 ) 


with 


r := Ut F(t) 

Which are, essentially the uncoupled equations of motion for the structure we 
are examining. 

1.1.2 Orthogonality Conditions 

By the transformation employed in eq. (1.1. 1-6) , we have used the so-called 
orthogonality conditions, that must be validated in order for the transformation 
to work. These are 


For the mass matrix : 


Uj> 

M u k = 0 

for all j ^ k 

(1.1. 2-1) 

Uj« 

M u j = mg en/ j 

for j=k 

(1.1.2-la) 


and for the stiffness matrix : 
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Ujt JC u.k = 0 for all j ^ k (1.1.2-2a) 

Ujt K Uj = kg en ,j for j=k (1.1.2-2b) 

For a mechanical explanation of this, we can consider that the mass-forces 
-o>;2 M uk of the k-th eigenfunction does not perform work on the 
displacements uj of another (different) eigenfunction j. The work being done on 

this eigenfunction is through the mass-forces of its own eigenfunction 
-C0j2 M Uj. The numerical value of m gen , that results from the front and back 

multiplication of the mass matrix by the modal matrix, is called the generalized 

mass . 


In addition, inserting (1.1. 2-1) and (1.1.2-2) into the equations for the eigenvalues 
we get an additional n equations, for the natural frequencies: 


-CDj 2 irigen,j + kgen,j ® 


(1.1 .2-3) 


or 

2 = ^gggj- (1.1. 2-4) 

) m gen j 

The uncoupled equations of motion take on the following form: 


m gen,j j + kgen,j Q j r j j 1 ,..., n (1.1 5) 

These n equations of motion for one DOF systems can be solved one after the 
other and after superimposing these solutions we get the total response of the 

system : 


x(t) = 



Uj qj = U q 


( 1 . 1 . 2 - 6 ) 


This equation also represents the transformation from the physical coordinate 
system to the modal coordinate system with the orthogonal eigenvectors of the 
structure as its basis. 
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A more intuitive and useful modal decomposition procedure is to choose the 
magnitude if the eigenvectors u in such a manner that the front and back 
multiplication of the mass matrix by the modal matrix results in a unity matrix 
for the generalized mass matrix and a stiffness matrix that has only co n 2 on its 

diagonal . 2 

<8> M © = I (1.1.2 -7) 


and 


© k © =ru 


( 1 . 1 . 2 - 8 ) 


The uncoupled equations of motion in modal coordinates now read 


I j + rx*J a = U T f 


(1.1. 2-9) 


With 

U T F =R 

as the generalized force and 
x(t) = © q 

as the transfer equation to modal coordinates. 


( 1 . 1 . 2 - 10 ) 


( 1 . 1 . 2 - 11 ) 


1.1.3 Structures with Distributed Mass 


2 Scaling the eigenvectors u to achieve this is termed orthonormalizing the eigenvectors. These 
orthonormalized eigenvectors are denoted <(> and found by 



The orthonormal modal matrix is ©. 
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In the case of a beam for instance, we cannot necessarily discretize the system as a 
lumped parameter model if we want to achieve a fair amount of accuracy. So, we 
take the mass distribution into account and set up a continuous system. The 
eigenvectors of such a system turn into eigenfunctions as the discretization 

becomes finer. 

The differential equations of motion for a beam can be quickly derived by 
applying the three kinds of equations used to solve any structural dynamics 
problem: equilibrium condition, material law and the kinematic law. Consider 
an element of the beam: 



H w dx 



dx 


Fig 1.1 .3-1 Beam Element 


Equilibrium Conditions 

Moment (disregarding 2nd order terms 

and the rotational inertia of the element) 

P)-o 

M” = Q 


Shear 


O' = p w -p 


1 0 


=> M" + p - |i W = 0 



Material Law 


M = El k (the linear relationship: bending moment- radius of curvature k) 


Kinematics 

k = -w" 

If the modulus of elasticity is constant, we can write 
B(x) = E I(x) 

and the differential equation of motion are 


[ B(x) w " (x) ]" + m(x) w(x) = p 
Or considering a uniform beam, the equation reduces to 

_ M //// <'**' +<* 

B w + m w = p 


1 

B(x),n(x) 

^ — — u. 

i 

^ \ » 

i Z,w(x) 


(1.1.3-1) 


(1.1.3-la) 


Fig 1.1. 3-2 Cantilever Beam 

As in the case of the lumped parameter model, we assume that it is also possible 
to modally decompose this beam into decoupled, one-DOF oscillatory systems. In 
addition we stated that the response of such a system could be expressed by a 
superposition of its eigenvectors, Uj, and a participation factor, q j, that carried the 

timevariation of the system. Now, in examining a continuous system, the 
eigenvectors are replaced by eigenfunctions £j(x). So the equation for the 

superposition of the eigenvectors now is 


w(x)= X^iM 3iM 


1 1 


(1.1. 3-2) 




In the case of a uniform beam it is still possible to find an analytical solution for 
the eigenfunctions. Using the boundary conditions for a cantilever beam, this 
solution is 


cosMXjx} - costXjx) sinh(X i x)-sin{A, i x} 
£i (x)= coshftql) + cos&il} ‘ sinh{X i l}+sinft. i l) 







’ 

" ^ 


First Eigenfunction 
or Mode 

X,l =1.875 


Second Eigenfunction 
or Mode 

jg = 4.694 


Third Eigenfunction 
or Mode 
\\ =7.885 


Fig 1.1. 3-3 First Three Eigenmodes 

Clearly, if the beam has a mass distribution and a changing geometry along its 
axis, an analytical solution can no longer be found and one must revert to 
numerics. A polynomial is used to approach the eigenfunction of the beam. The 
beam is discretized into small parts in a consistent manner and eigenvectors of 
the discrete system are calculated. The element of the eigenvector u n is 
understood to be the deflections of the structure at that lumped mass m nn The 
polynomial is laid through these points of deflections in order to attain a 
estimate of the shape- or eigenfunction. 


The polynomial approach is written 

*W»S c ij Tlj( x ) 
j 

inserting this into eq. (1.1. 3-2) gives us: 


(1.1. 3-4) 


w(x)= X X c ij HjM 3iW (1.1. 3-5) 

' i 

where qj being the coefficient matrix for the polynomial approach. In closing, 
note the similarities between the discrete formulation and the continous 
formulation of the problem: 


m j qj +d j qj +k j qj = 


(1.1. 3-5) 






Discrete System 
i= 1...N 

Discrete System 
j= 1...00 

Eigenforms or 
mode shapes 

Uj - vector 

£j(x) - function 

generalized mass mj 

ujt M uj 

| |i(x) (pj 2 (x) dx 
o 

generalized damping dj 

km m i +k s s. 

kminj+k.. Sj 

generalized stiffness kj 

uj* K uj 

/ B(x)^j" 2 (x)dx 
o 

generalized forces Pj 

uj* p f(t) 

j £j(x) p(x) dx 
o 

Response of System 

x = u ; q j 

j=i 

icF 

>< 

sWi 

II 


Table 1.13-1 Continous and Discrete Systems 


1.1.4 The Finite Element Theory 

The Finite Element Theory is an approach that discretisises the system into 
small, well defined pieces whose properties are known, and derives the vibration 
characteristics by assembling these small parts back to the large system and 
solving the equation of motion and decoupling them. The use of finite element 
techniques for dynamic structure analysis introduces great flexibility to rotorcraft 
design which is needed because of the often complex root geometry and 
nonlinear stiffness and mass distribution of the blade. The application and 
implementaion of Finite Elements will be discussed in detail later. 


1.2 Experimental Modal Analysis 




Experimental measurements are usually made in the Laplace-Domain (also 
known as the Frequency Domain ) . This form of vibrational analysis became 
popular as very fast algorithms for the Fast Fourier Transform (FFT) and the 
Discrete Fourier Transform (DFT) were developed in the 1960's. The ability of 
these Digital Fourier Analyzers to quickly convert time histories and extract 
modal characteristics in the frequency domain has made them an important tool 
in measuring the modal characteristics of structures. The NASA shake test 
identified the modal residues of the UH-60A rotor blade. To correctly interpret 
this data we must look a little closer at the theory behind modal residues, the 
FFT and the Laplace-domain. 

Consider the second order differential equation of motion derived in a previous 
chapter : 

M x +D x +K x = F(t) (l.l.l-l) 

1.2.1 The Laplace Transform 

Any function of time can be transformed into an analytical function of the 
complex variable, s, in the Laplace domain by the Laplace transform. There it can 
be algebraically altered and transformed back into the time domain. The greatest 
advantage that this transformation gives is that the differential equation of 
motion is changed into a simple polynomial that is much easier to solve. The 
transformation equation to the Laplace-domain from the time-domain is 

oo 

J C [x(t)] = X(s) = | e-stjxft)] dt (121-1) 

o 

and the back-transform is 

OO 

L _1 [X(s)] = x(t) = J e* st [X(s)] dt (1.2. 1-2) 

o 

The complex Laplace variable s gives the relationship and magnitude of 
frequency and damping parameters on a complex plane is shown in Fig 12.1-1. 
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I,a place Domain 



The transform of a variable x to the Laplace domain will be written as JL lx] . To 
use this Laplace transform on our differential equation we need the 

transformation of x and x . 


4f]- H#* 


Integrating this by parts finally gives us 
X[f] = sX[x]-x(0) 


(1.2. 1-4) 


L [x]=s X(s)-x(0) (1.2.1-5) 

Similary we can find 

£,[§£[• L [ x] = S 2 X(s) - s x(0) - x(0) (1.2.1-5) 

If we ignore the initial conditions (which can always be done for stable systems) 
this transforms our equations of motion to 




M s 2 X(s) + D s X(s) + K X(s) = F(s) 


( 1 . 2 . 1 - 5 ) 


an algebraic equation in the complex variables. It can be easily seen that the 
mathematics do indeed take on a simpler form in this domain. 

Solving for X(s) 


X(s) = 


F(s) 

B(s) 


( 1 . 2 . 1 - 5 ) 


where X(s) = Laplace Transform of the Responce 

F(s) = Laplace Transform of the Excitation 

B(s) = M s 2 + D s + K 

This equation is used as the basis to define the Frequency Response Function 
(FRF) such that: 



L.-transform of the Output = FRF * 

L.-transform of the Input 

or 

X(s) = H(s) 

F(s) 


with H(s) = [ M s 2 + Ds+ K ] 1 , the Frequency Response Matrix 

Now all the outside forces, responses of the structure, velocities, accelerations 
and coordinates are expressed in the frequency-domain. The transfer matrix H 
contains transfer functions h(s) that give the transformed response X(s) for every 


input f(s) at every DOF. 



Xj(s) 

X 2 (s) 


hn(s) hi 2 (s) • • • h-j n (s) 

fl(s) 

f 2 (s) 


X 3 (s) 

— 

h 21 (s) 

f 3 (s) 


_ X n (s) _ 


* * • • • ♦ 

_ h n i(s) • • • • - 

_ W _ 


1 6 


The transfer functions, h(x), are also of a complex nature. The real part 
corresponds to the magnitude and the imaginary part corresponds to the phase 




shift of the response. For simplicity, let us once again consider a simple two-DOF 
system, keeping in mind that the statements made can be transferred to 
multiple-DOF systems as well. 

In this case, the characteristic equation (CE) for the inverse of the frequency 
response matrix B would be a polynomial of the degree 2n with n pairs of 
conjugate complex roots, p k and p k * . Expressing the CE in terms of its roots we 

can write: 

Det B = C (s- Pl ) (s-pi*) (s-p 2 ) (s-p 2 *) •••• (s-pn> (s-p n *> (1.2.1-6) 

So, fl\e components of the frequency response matrix can be written as 


hjj(s) 


m,j S 2 + djj s + k jj 
CE 


(1. 2.1-7) 


A diagram representing the properties of the FRF and its connections to the 
DOFs is shown in Fig. 1. 2.1-2. The input fj as wll as the outputs x-j ... x n are 
transformed (=>) into the Laplace domain becoming Fi and Xi ... X n . The 
magnitude and phase shift of the output in relationship to the input excitation is 
given by the FRF. 



Fig 1. 2.1-2 The Frequency Response Function 



If we assume that the poles of H, i.e. the roots of the CE of B, are of multiplicity 

one for our physical system, we can expand the denominator 3 and write 
n 

h ij (s)= [(S-Pk) + (s-Pk*) (1 ‘ 21 8) 

k=l 

Pk=Ok + icok 

Ok = damping coefficient 

C0k = nat. freq of oscillation 


Examining the responce matrix for one specific mode, k, the response functions 
hjj all have the common denominator |^ s _p k ) + (s-pk*)J differ only in the 

numerator, rk . These residues rk are the only parameters that vary along the 
structure and are collected in the residue matrix Rk* 


i'll r 12 * * * r ln 


L rnl J 

implying that H = or H k ~ R k 

The matrix of the system residues, just like the modal matrix, represent the 
motion of a structure when excited at a resonant frequency, both giving the 
amplitude and direction of the response at some DOF. Taking advantage of the 
symmetry inherent in the response matrix H (hq = hjj) we can make the 

connection to the modal matrix. This represents the experimental /analytical link 
of the modal analysis procedure. 

R k = Qk Uk Uk T (1.2.1-10) 


3 A polynomial B(x) that is divided by a polynomial A(x) can always be written in the partial 
fraction form 


A(x) _ r l 
B(x) " 3^ 


+ 


S -Pn 


+ k(x) 


(s-pi) being the roots of B 


r 2 

*P2 



where is an arbitrary constant since the mode shapes are defined only in 
direction not magnitude. 


Arbitrarily choosing the i-th column of the residue matrix, we can write 
(supressing the index k for the moment, since we are examining only the 
transfer function for the k-th mode) 


— — 


— 

r li 


Ui Ui 

*2i 


u 2 Ui 

r 3i 


u 3 u i 

rii 

= Q 

Uj Uj 

_ r ni _ 


_ u n u i _ 


= Q Uj u 


( 1 . 2 . 1 - 10 ) 


We can calculate the residue for the driving point (structure excitation point), r H 


r ii = Q u i u i =QUj 2 


( 1 . 2 . 1 - 11 ) 


or 

( 1 . 2 . 1 - 12 ) 

which delivers the proportionality factor of the two matrices. Thus, as long as a 
driving point measurement is made and this driving point does not correspond 
to a node of the examined mode, the entire residue matrix can be constructed by 
measuring only one row or column of the transfer matrix, i.e. knowledge of one 
mode. 




2. Modeling the Beam by the Finite Element Method 
2.1 Introduction 

Since a wealth of data concerning the vibrations of a cantilever beam exists 
(Ref. 2.-1,6,11,12), a finite element model for the boundary conditions clamped- 
free was written. Once a sufficient accuracy was determined, this model could be 
altered to encompass other boundary conditions of a vibrating system, i.e. 
hinged-free or free-free. 

Because the algorithm for predicting the rotating mode shapes (see section 
2.3.1) relies heavily on the slopes or derivatives of the eigenform of the non- 
rotating mode shape of the beam, these slopes were compared with those given 
in Ref. 2.-1 and Ref. 2.-2. After verifying the accuracy of the non-rotating model, 
an extension of the method of Ref. 2.-3 using the finite element method was 
applied, and the results for the rotating model were compared to the exact 
solution of Ref. 2.-2 - a 5 mode extension for a cantilevered Beam and a 6 mode 

extension in the case of a hinged beam. 

The beam is divided into n elements and a stiffness and mass matrix is 
derived for each single element. These element matrices are then assembled into 
a system matrix, taking into account the boundary conditions. In this case the 
continuous mass distribution is replaced by a lumped parameter mass 
distribution, which means that each end (or node) of the element is assigned half 
the finite element mass. 

2.1.1 The Element Stiffness Matrix 

The element stiffness matrix, k, expresses the relationship between static 
forces and moments, P, and linear and angular displacements, & • 

P=k 5 (2.1.1 1) 

A beam segment or element with the properties shown in Fig. 2.1.1— 1 has four 
degrees of freedom, a linear and an angular one at each node. 
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Nodes 



Dftnrees o f Freedom 
u linear displacement 
<P angular displacement 


Fig. 2. 1.1-1 BEAM Element 


The element stiffness matrix is derived in detail in Ref. 2.-4 using the principle 
of virtual deflection and found to be 



12 -61 -12 -61 
-61 4l 2 61 l 2 
-12 61 12 61 
_ -61 l 2 61 41 2 _ 


( 2 . 1 . 1 - 2 ) 


(4x4) 

2.1.2 The Element Mass Matrix 

Following the lumped mass distribution assumption, the inertia 
associated with each rotational degree of freedom is assumed to be zero while the 
mass of each element is assigned evenly to the nodes. This distribution of the 

mass is shown in Fig. 2.1 .2-1. 

m 

2 2 
Fig. 2.1. 2-1 Lumped Mass Distribution 

Thus, the element mass matrix can be written as 

y 0 00 

0 0 0 0 
m = m 

0 j 0 0 
_ 0 0 0 0 _ 
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2.1.3 Assembly of the Structure Stiffness Matrix 

The stiffness matrix of the structure, in our case a beam consisting of n 
finite elements, is assembled by identifying the degrees of freedom at each 
element node, numbering them in a consistent manner and simply adding those 
contributing to the same system degree of freedom. This is what is known as the 
direct method. For an homogeneous cantilever beam with a constant cross-section 
(we will call this sort of beam simply "constant"), this assembly process is 
demonstrated in Fig. 2. 1.3-1. For a cantilever beam, the boundary conditions are 
<pi = 0 and uj = 0 . 



Fig. 2. 1.3-1 Assembly of the Stiffness Matrix 


The system stiffness matrix relates the forces and the displacements at the system 
nodal coordinates in the same way as the element stiffness matrix relates them in 
the element nodal coordinates. The rules governing the assembly process for the 
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system matrix are geometric compatibility at the element nodes, i.e. the 
displacement at the nodes shared by the elements must be the same for each such 

element. 

2.1.4 Assembly of the System Mass Matrix 

The assembly of the mass matrix is a simple matter of adding the 
contributions of the lumped masses at the nodal coordinates. Since no inertia is 
assigned to the rotational degrees of freedom, no mass terms are found in the 
rows and columns associated with these degrees of freedom. Following this, the 
positions of the mass matrix referring to the rotational degrees of freedom are set 
to zero. This assembly procedure is shown in Fig. 2.1. 4-1 for a beam consisting of 
10 finite elements and thus has 20 nodes. 



Fig. 2. 1.4-1 Assembly of the Mass Matrix 

2.1.5 Static Reduction of the System Matrices 

Subdividing a structure into many finite elements leads to very large 
matrices. In the case of a simple beam divided into n segments, the size is [ 2n x 
2n 1 . A way of reducing the size of the structure matrix is to identify the degrees 
of freedom that are not needed or those that are not of interest and ascertain 
their dependency upon the remaining degrees of freedom. This is called static 
condensation. In the case of lumped masses, where no inertia is assigned to the 
rotational degrees of freedom, the following condensation technique not only 
leads to a reduced mass and stiffness matrix, but also leads to an equivalent 
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eigenproblem. The reduction of massless or inertialess degrees of freedom 
introduces no error. 

The reduction of the stiffness matrix is somewhat involved, but can 
nonetheless be easily accomplished. First, we permute the rows and columns in 
the stiffness matrix such that the displacement and rotational degrees of freedom 
are adjacent to each other. This leads us to a stiffness matrix with rearranged 
rows and columns, K re that can be partitioned as follows: 

ul 
u2 


♦1 

* 

' (2.1. 5-1) 


K = 

re — 


K 


uu 


K 


up 


K 


I K 


<pu 


<W> 


So that now the equations of motion can be written as 


r~ — 


P — 


— — 


•— — 


— — 

m o 


•• 

U 


K K 


u 


p u 





UU Uf 





0 0 


•• 

<p 

+ 

K K 


<p 


0 

_ — 


— — 


— — 


— — 


— — 


The lower part of this matrix equation gives us 
K 9U u +K w £ =0 


(2.1 .5-2) 


(2.1. 5-3) 


[— 




'o' 

UJ 


K w 


.0. 


(2.1. 5-4) 


as the equation of constraint. Thus, the rotational and translational degrees of 
freedom can be related to each other as 

(J> = -K w -1K 9U u 
Setting 

Tcpu = " ^<p(p ^ ^<pu 

the transformation equation becomes 


(2.1. 5-5) 
(2.1. 5-6) 
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1 

L«J 

Now, the rearranged K re and M re matrices can be reduced in dimension by the 
transformation matrix T to 


L <P U J 


u 


= T u 


(2.1.5-7) 


M =T l M re T 
K =T‘ K re T 


(2.1.5-8a) 

(2.1.5-8b) 


The mass matrix is reduced, in the case of lumped masses, simply by 
deleting those rows and columns that pertain to the rotating degrees of freedom. 


or 


A re 


K uu — " 


IU 


-1 


(2.1.5-9a) 
(2.1.5-9b) 

The rearranged mass matrix, M re , remains essentially unchanged by the 
reduction, while the stiffness matrix is transformed to the reduced stiffness matrix 
K uu . The eigenvalue problem 


A 

(K 


k UU 


re 


) & = 0 


(2.1.5-10) 


can now be solved using the available algorithms. 


2.1.6 Treatment of Singular Matrices 

A beam with one end free and the other end hinged represents a statically 
undetermined system, and the system stiffness matrix assembled using the finite 
element method is singular. A singular matrix cannot be inverted. In addition, 
the numerical size of the elements of the mass and stiffness matrix differ by a 
large factor that increases with the degree of discretisation. The above factors 
contribute to the danger of numerical difficulties in solving the eigenproblem. 
To circumvent these numerical problems, the mass and stiffness matrices were 
divided by their norm, and the equivalent eigenproblem was solved to obtain 
the eigenvectors and eigenvalues . 
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The norm of a matrix A is defined as 

DA II = (det ( A At ))V2 


( 2 . 1 . 6 - 1 ) 


Dividing a matrix by its norm in this case yields a numerically well-conditioned 
matrix. Applying this to the eigenproblem results in 

det [ K' - T| M' ] = 0 (2.1. 6-2) 


with 



M . II Mil 

M ' = 11 M II 11 “ * UK II 


(2.1. 6-3) 


After solving for X, we reinsert the eigenvalues (one of them being zero; 
corresponding in this case to the rigid body mode of the hinged beam) into the 
eigenproblem (2.1.5-10). The n-th eigenvector can now be found by inserting the 
n-th eigenvalue into the eigenproblem: 


Inserting = X n into 

( K- X[ I M ) <J»i = A 4>| =0 (2.1 .6-4) 

where <j>j is the eigenvector corresponding to X[ t and decomposing A into the 
matrices P and P ^ and the diagonal matrix D 


A = P D P - 1 


(2.1. 6-5) 


we can write 

DP' 1 4> = 0 

Defining 

y_ := P'l <J> 

it follows that 
Dy_ = 0 


(2.1 .6-6) 


(2.1 .6-7) 


(2.1 .6-8) 
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D is a diagonal matrix containing the eigenvalues of A. The matrix A is of 

course singular, since an eigenvalue was inserted, so it follows that one or more 
of the diagonal elements of D are zero. The elements of y that are in the same 
rows that in D contain zero can be arbitrarily selected in order to satisfy the above 
equation. Correspondingly, all elements of y in the same rows that in D are non- 
zero must contain a zero in order to satisfy the equation. 

Thus, the following equation can be constructed: 


D _y 






0 


with yi— ym-i arbitrary 

and y m -y n =0 


( 2 . 1 . 6 - 9 ) 


So the eigenvectors <J>i , corresponding to the eigenvalues Jlj , can be determined 
by simply solving (2.1. 6-7) in the form 

* =P * ( 2 . 1 . 6 - 9 ) 
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2.2 Results of the Non-Rotating (Static) Analysis 


Calculations were performed on cantilever and hinged beams, having 
both constant and tapering stiffness and height distributions. Results are 
presented for tapered and constant cantilevered and hinged beams in Figs. 2.2-1 
to 2.2-4. 
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Fig. 2.2-1 Mode shapes and Coordinates of the Cantilevered, Constant Beam 
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Fig. 2.2-2 Mode shapes and Coordinates of the Cantilevered/Tapered Beam 
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Fig. 2.2-3 Mode shapes and Coordinates of the Hinged, Constant Beam 
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Fig. 2.2-4 Mode shapes and Coordinates of the Hinged, Tapered Beam 
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2.3 Prediction of the Rotating Natural Frequencies 

2.3.1 Predicting the Bending Frequencies of a Rotating Beam given 
the Non-Rotating Bending Frequencies. 

To verify the accuracy of the mathematical rotor blade model described in 
section 2.2. , we will compare the results of that simulation to those obtained by a 
different approach for beams with well known geometric properties and later for 
the Blackhawk rotorblade itself. Since the non-rotating bending frequencies of 
the Blackhawk rotor blade with free-free boundary conditions are known from 
the shake test (Ref. 2.-7), we will use a Lagrangian method suggested by Ref. 2.-3 
to predict the rotating frequencies and modes. This approach allows us to 
determine the rotating bending frequencies for any rotational speed, given only 
the structural properties of the rotor blade or beam that is to be analyzed. There is 
no restriction as to the stiffness or mass distribution. Only the stiffness and mass 
matrices must be assembled in order to determine the non-rotating mode shapes 
or eigenvectors. Both Loewy (Ref. 2.-3 ) and Yntema (Ref. 2.-2) make the 
simplifying assumption that a linear height distribution leads to both a linear 
mass and a linear stiffness distribution. Using the finite element approach will 
allow us to negate this simplification and calculate beams or systems with a truly 
linear height distribution ( resulting in a cubic stiffness distribution ) and systems 
with a truly linear stiffness distribution (resulting in a third-root height 
distribution). As we will see, the calculations using this more exact approach will 
im prove upon the accuracy of the previous results. 

The bending deflections of the rotating mode will be expressed in terms of 
the deflections of the non-rotating modes 

Wj(t) = X Ui (n) q(t) n (2.3.1-1) 

n 

where 

Wj(t) = the i-th rotating mode shape vector 

Ui<"> = the n-th non-rotating mode shape 
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q(t) n =the time and amplitude dependency of the n-th non- 
rotating mode (and thus also the n-th generalized 
coordinate in the Lagrangian approach ) 


Following (Ref. 2.-3) , the kinetic energy terms are expressed as 

T = | * m i) (2.3.1-2) 

and after the differentiations demanded by Lagrange's equation 


a 

at 


fdT } 




(2.3.1-3) 


these kinetic terms turn out to be 


Mgen & 


(2.3.1-4) 


where Mg en is the generalized mass matrix obtained from the static (non- 
rotating) analysis of the beam, and q 4 is the second derivative with respect to 
time of the generalized coordinate. 

The potential energy terms are divided into two parts. The first, Uj, 

represents the potential energy stored in the beam due to deformation (elastic 
properties) and the second, U 2 , represents the work done by centrifugal forces on 
the lumped masses, m,, acting through changes in the radial position, xj of these 
masses as the blade bends. This "shortening" A* of the blade at the radial position 
xj along the beam coordinate x 0 can be expressed as 


Ai = 



(2.3. 1-5) 


The expression for Aj is derived and discussed at length in Ref. 2.-5. The 
shortening can be approximated with the the following summation: 


33 


(2.3.1-6) 



where j = 1 is the position at the root of the beam. 

Inserting the above expressions into the potential energy terms 


U! - 2 5-r ®^gen 3^ 
n 

{ 23 . 1 - 7 ) 

U 2 = Zm i x oi a 2 Aj 
i 

{ 23 . 1 - 8 ) 

and considering that 


Kg en = 0)2M gen = S m i w i (n)2 

{ 23 . 1 - 9)1 

we obtain the equation for the r-th generalized coordinate (the r-th row of the 

rotating beam equation) as 


~ V ( ^ w(n) '\ 2 ( 3w<r) ' 

Mgen (r) % + x oi / ,qr| 3^ 3^ 

j =1 

v 2 

1 AXqj + ^gen qr ~° 

(2.3.1-10) 

Substituting 

i 

Tm - [ 3*o 1 ^ 

1 H 

yields 

(2.3.1-11) 

Mg en q + n 2 Yq + K gen a = 0 

(2.3.1-12) 


The term K stat = K gen is the static or non-rotating part of the beams 
stiffness matrix, and K dyn = 0 2 Y is the dynamic part of the stiffness matrix 
brought on by the centrifugal forces on the lumped masses. 

Adding K stat and K dyn leads to the well-known problem 


1 if the condition that M is diagonal holds 


M <j + K 3 = 0 


(2.3.1-13) 


that can be solved by the usual algorythms to obtain the rotating beam modes. 

Taking into account the first 5 non-rotating beam modes of a beam subdivided 
into 20 lumped masses, we obtain a K dyn that is populated in the upper left hand 

corner (5x5) with the elements 

Yij (with i,j=1...5) 

2.3.2 Linear Height and Stiffness Distributions 

Calculations to verify the algorithm prior to applying it to the UH-60 rotor 
blade were done for the following examples in order to gain some insight into 
the order of accuracy. These results were compared to the results of a 5-mode 
extension performed in Ref. 2.-2 for cantilever and hinged beams of constant and 
linear varying stiffness distributions. As stated earlier, the five and six mode 
extensions were based upon the assumption that a linear height distribution 
implied a linear mass and a linear stiffness distribution. This analysis introduces 
results for truly linear distributions of the above parameters. For example, the 
results obtained in Ref. 2.-3 are for a beam of linear varying height taken from the 
book by Bisplinghoff Ref. 2.-6 whereas the results are compared to those of 
Yntema (Ref. 2.-2), who only examines beams of a constant or linearly varying 
stiffness distribution. The height function of the beam found in Ref. 2.-6 is 

h(x)=2 b(l-X*) (2.3.2-1) 

with A.=0.2 . Therefore the stiffness is : 

1 = h b h3= h. b2 (1 ' X 3 (2.3. 2-2) 

X 

where x = ^ is the dimensionless blade coordinate. 

The stiffness will only be linear if the higher terms of x can be dropped. 
This cannot be supported for a value of X = 0.2 . Thus, a height-variation 
producing a truly linear stiffness distribution was used in our examples. Also, in 
the hopes of obtaining a better estimate for the rotating modes, the first five non- 


35 


rotating mode shapes were included, instead of only the first two as in Ref. 2.-3. 
Because the algorythm uses the slopes of the non-rotating mode shapes, the 
possibility of using a polynomial to obtain the derivatives instead of numerically 
calculating the slopes by a three-point and five-point differentiation rule was 
addressed, since a sufficiently exact polynomial can easily be fitted to the first five 
mode shape vectors. This did, in fact give a better approximation of the slope of 
the beam at the lumped mass in question, but, as pointed out in Ref. 2.-3, the 
number of mass stations makes a far greater contribution to the convergence of 
the algorythm towards the correct solution than improving the accuracy of the 
derivatives. 

2.3.3 Determining the Rotating Mode Shapes 

After having determined the rotating frequencies by the addition of the 
dynamic stiffness matrix expressed in terms of the generalized system equations 
it becomes necessary to do a back-transformation prior to obtaining the mode 
shapes for the rotating system in the usual Cartesian system coordinates. We will 
use lowercase to indicate the generalized coordinates and uppercase to indicate 
the system (Cartesian) coordinates. 

We can express the rotating system stiffness matrix, k S y St/ in generalized 
coordinates 

Kyst = kgen + k dyn (2.3.3-1) 

The diagonalizsed stiffness matrix,k gen , was found by transforming the static 
stiffness matrix, K stat , by the system modal matrix, U, a matrix composed of the 
eigenvectors 2 obtained while solving (2.3.1-13) 

kgen =U T K stat U (2.3 .3-2) 

Inserting (2.3.3-1) into (2.3.3-2) we get 


2 since U is composed of the eigenvectors of our vibratory system arranged in columns, and 
these eigenvectors have been normalized prior to placing them into U, we have an 
orthonormal matrix. The transformation equation U' 1 IK ID can be simplified to IK U, 
since U’ 1 = U T for orthonormal Matrices. 


36 


(2.33-3) 


^syst — ^ ^stat ^ + ^dyn 
Transforming (23.3-3) into Cartesian coordinates with the system matrix, U, 
gives us 

(uV k syst (U)' 1 = K stat + (uV k dyn (U)- 1 (23.3-4) 

or 

K syst = Kstat + K dy „ (23.3-5) 

To find the rotating eigenvectors, we must solve the equation 


M q + K syst q = 0 


(23.3-6) 


that leads to a new system matrix, U rot , that is composed of the rotating 
eigenvectors. This new system matrix diagonalizes the rotating stiffness matrix. 


^rot - kgen,rot 


(2.33-7) 


The rotating and non-rotating normalized deflection curves of a 
cantilever constant beam for the l s ^ to 3 r< ^ natural frequencies are shown in Fig. 
23.3-1. For the hinged constant beam these are shown in Fig. 2.33-2. The 
rotational parameter 

(2.33-8) 

was set at 4 for the cantilever beam and at 2 for the hinged beam. It is seen that 
the discrepancies between the non-rotating and the rotating mode shapes are 
small, and that the nodes are virtually coincident. 
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2.3.4 Natural Frequency as a 


Function of Rotational Speed 


Figs. 2.3. 4-1 to 2.3.4-4 give an overview in how the natural frequencies 
increase with increasing rotational speeds. The factor ¥ relates the square of the 
first non-rotating frequency C0 nl to the square of the frequency that the system is 

rotated at, Q . . This parameter stays the same, regardless of which mode is being 


examined. 




(23.3-7) 


The factor for the abscissa,C>, shows the magnitude of the change in 
rotating to non-rotating natural frequency. The denominator of this parameter 
always holds the non-rotating natural frequency of the mode being examined. 


w 


ru 


= ^rot i¥ (2.3.4-1) 

~l<°nij 

The parameters are squared in order to better asses the reliability of the obtained 
results when comparing them to Ref. 2.-3. Notice that squaring the parameter 
also squares the error. The diagram for the cantilever beam with non-constant 
height distributions (Fig. 2.3.4-2), reaffirms the previously stated concerns. Only 
when a beam is considered that has a linear height distribution and a linear 
stiffness distribution are examined, does one receive the results stated by Yntema 
(Ref. 2.-2). This form of beam, having both a linear mass distribution as well as a 
linear stiffness distribution does not exist for isotropic structures. 

The frequency graph for the "tapered" cantilever beam (Fig. 2.3.4-4) shows 
the frequency increases for a beam of actual linearly decreasing height and 
linearly decreasing stiffness at a constant width. 
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2.3.5 Accuracy and Convergence Characteristics for this Finite 
Element Application 

It is of interest to us how exact our eigenvalues and eigenvectors are, and 
how confident we can feel about our results. The solutions we have found for 
the calculated ( ~ ) eigenvectors, $ , and for the calculated ( ) eigenvalues, X , 
were calculated by the EISPACK matrix and eigensystem routines (Ref. 2.-8). Of 
interest would be the convergence of the finite element program as a function of 
the number of elements used in the discretisation of the beam. 

Having made the rather bold assumption that the rotational degrees of 
freedom were not assigned any rotary inertia, it is of even bigger importance to 
know more about what effect this simplification has on the calculations m order 
to be able to judge the effect on the projected rotating modes and frequencies. 

Throughout literature ( e.g. Ref. 2.-1, Ref. 2.-4, ...) analytic results are found 
for the mode shapes for constant beams with a constant stiffness along the beam 
under the assumption of a variety of boundary conditions. We will compare our 
solutions with the analytic solutions for the cantilever and hinged beams and 
attempt to show the convergence characteristics for the first three non-rigid body 
eigenvalues as a function number of elements used. A strategy for determining 
the quality of a calculated eigenvalue even if the true eigenvalue is unknown is 

also introduced. 

Exact analytic solutions for the hinged-free, free-free and cantilevered beam 
natural frequencies coj^ are shown in Table 2.3. 5-1 using 

a3 ' M> 


Table 2.3.5-1 Exact Solutions for the Constant Beam 


System Type Eigenvalue Equation Rigid Body Values for Xk 1 

Modes 


Cantilevered 1+cos X\ ch X\ = 0 
Hinged-Free (XI) (th Xl-tan XI) = 0 
Free-Free (XI) 4 (1-chXl cosXl) = 0 



1,8751 / 4,69409 / ... / (n-0.50)rc 
3,9266 / 7,06858 / ... / (n+0,25)rc 
4,73004 / 7.8532 /... / (n+0,50)7i 
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The convergence characteristics of the implemented finite element algorithm 
toward the exact analytical solution are shown in the following figures. Figs. 

2.3.5- la to 2.3.5-c show the convergence of the calculated solution divided by the 
exact solution of the cantilever beam eigenvalues towards the normalized exact 
solution (1.0) as a function of number of elements used, while Figs. 2.3.5-2a to 

2.3.5- 2d show the same for the hinged beam. 


Cantilevered Beam 



Fig 2.3.5-la 1st Natural Frequency 



Hinged Beam 


The 1st natural mode 
for the hinged beam is 
the rigid body mode 


Fig 2.3.5-2a 1st Natural Frequency 
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Pig 2.3.5-lc 3rd Natural Frequency Fig 2.3.5-2C 3rd Natural Frequency 



Number of Elements 


Fig 2.3.5-2d 4th Natural Frequency 
The analytical equation for the mode shapes is 


ch X-jx - cos XjX sh XjX - sin X.jX 
w i (x) = ch X.jl + cos X-il ‘ sh Xil + sin X.jl 


(23.5-2) 


for the cantilever beam and 

sh X.jl sin Xjx + sh X.jX sin X.,1 
w i<*) = 2 sh Xjl sin X? 


for the hinged beam 

Figs 2.3.5-la to -2d show that a good accuracy can be obtained by using 20 
finite elements on a beam, resulting in having to solve only a 10x10 matrix after 
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static condensation. Tables for the first few eigenvectors for the calculated beams 
can be found in Appendix A. 

2.4 Comparison of the Finite Element Calculation with the Full- 
Scale UH-60A Rotor Blade Non-Rotating Modal Analysis Shake 

Test 


2.4.1 NASA Shake Test 

A shake test of the UH-60A rotor blade was performed at the NASA Ames 
Research Center during April 1986 as part of the Phase I Flight Test MTR 
(Modem Technology Rotor) project. The rotor blade was outfitted with strain 
gages along the span and a movable piezoelectric accelerometer was used to 
measure the response at the node points. The data was analyzed to obtain the 
frequency response function (FRF) of the blade excitation and response. The 
linearity of the rotor blade system was established by varying that the FRF 
remained constant for an increasing input force level. The minimal input force 
level for a constant FRF was found to be 1 pound-force ( 4.44822 N ). An upper 
boundary can be set at a force level that produces large, non-linear deflections. 
The Maxwell-Betty theorem was also found to apply. This states that the FRF that 
gives a response at point i of the structure resulting from a force input at point j 
is the same one as for a response at point j resulting from an input at point i (Hij 
= Hji) as shown in Fig. 2.4.1-1. This relationship was tested at frequencies of 12.55 

Hz and 25.11 Hz and shown to hold . 


Input Exitation : e 
System Response : r 



Fig. 2.4.1-1 Maxwell-Betty Theorem 
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The rigid body modes for the rotor blade suspended in a pseudo free-free 
manner by bungee chords was found to be 0.8 Hz for the 1 st flapping and 0.2 Hz 
for the 1 st torsional modes. 


Trailing Edge 


35 33 
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P 





□ 
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□ 


Leading Edge 

Fig. 2.4.1-2 Accelerometer Measurement Positions 


The FRF and the modal data were measured by placing the accelerometer 
at coordinates 16 in (40.6 cm) apart along the leading and trailing edges and at the 
node points of each identified mode. The accelerometer measurement positions 
are shown in Fig. 2.4.1-2. The blade was also examined by CAMRAD 
(Comprehensive Analytical Model of Rotorcraft Aerodynamics and Dynamics) 
in a hinged-free configuration and these results were transfered to the shake test 
free-free boundary conditions by calculating the ratios between the exact free-free 
and hinged-free solutions and applying them to the CAMRAD hinged-free 
results (Ref. 2.-7). This may account for some of the large errors found with the 
CAMRAD program, especially in the higher modes, because the ratios were 
found by looking at the vibrating beam partial differential equation 


dx 


EI(x> 


d 2 y(x,t) 

dx 2 


= -m(x) 


9 2 y(*,0 

2 

dt 


(2.4.1-1) 


for the two boundary conditions and then calculating the ratios. Since the above 
equation applies only to beams of a constant stiffness and mass distribution there 
is some error introduced by using it to obtain the ratios of hinged-free to free- 
free. The shake test identified 5 flapwise modes along with 2 chordwise and 2 
torsional modes. The 5 flapwise modes that are examined and compared to our 
finite element model are given in the form of modal residues at the leading and 
trailing edge. The residues (given in Appendix B) represent the amplitude of the 
mode at the point of measurement. The residues of the trailing edge differ from 
those of the leading edge since the individual modes are not totally uncoupled. 
To negate this effect we will calculate the residues at the quarter-chord. 
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supposing that the residues change in a linear fashion across the span of the 
blade as demonstrated in Fig. 2.4.1-3. When these residues are normalized they 



SC1095-R8 

Airfoil 

Fig. 2.4.1-3 Residue Calculation 
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2.4.2 Shake Test Results 

2.4.2.1 Rotor Blade Data for Shake Test 

The modulus of elasticity, E, for the Blackhawk rotor blade is 
E = iq 6 [ lb 2 -j = 0.68 10 10 The stiffness and mass distribution are given 

in Table 2.4.2.1-1 and in Figs. 2.4.2.1-1 and 2.4.2.1-2 and the natural frequencies 
for out of plane flapping that were measured during the shake test as well as the 
values calculated by CAMRAD are given in Table 2.4.2.1-2. 


Table 2.4.2.1-1 Mass and Stiffness Distributions of the Blackhawk Rotor Blade 
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1.00 


Root Blade Coordinate [cm] 

Fig. 2.4.2.1-1 Stiffness Distribution of the Blackhawk Rotor Blade 


Root Blade Coordinate [cm] 1 

Fig. 2.4.2.1-2 Mass Distribution of the Blackhawk Rotor Blade 
Table 2.4.2.2-1 Natural Frequencies, Out-of-Plane Flapping 


Mode 

Shake Test [Hz] 

CAMRAD [Hz] 

1st 

4.34 

4.67 

2nd 

12.55 

12.75 

3rd 

24.99 

25.49 

4th 

41.63 

42.32 

5th 

63.71 

57.75 






2.4.3 Results of the Non-Rotating Finite Element Analysis 

2.4.3.1 Discretization Technique 

The Blade was divided into three large parts in accordance with the degree 
that the mass and stiffness distributions varied. These parts were then in turn 
subdivided into finite elements. Three different division schemes (DS) were 
introduced to investigate the way the solutions converge as a function of tire 
number of elements. The DS are given in Table 2.4.3.1-1 and DS2 for the first 
normalized cantilever mode is shown in Fig. 2.4.3.1-1. 


Table 2.4.3.1-1 Division Schemes 


A UV1V *• 

1 

Scheme 1 

Scheme 2 

Scheme 3 

Number of Elements 
Parti (0...2m) 

20 

25 

35 

Number of Elements 
Part 2 (2.. .6m) 

15 

20 

30 

Number of Elements 
Part 3(6...7.29m) 

20 

25 

35 
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Results for the different DS and their deviations from the shake test results and 
CAMRAD are given in Table 2.4.3.1-2 below. 


Table 2.4.3.1-2 Results for Different DS 


Mode 

Nat. 

Frequency 
(Shake Test) 
[Hz] 

l auic j. 

Nat. 

Frequency 

(DS1) 

[Hz] 

Nat. 

Frequency 

(DS2) 

[Hz] 

Nat. 

Frequency 

(DS3) 

[Hz] 

Change: 

DS3(FEM)- 

Test 

Change: 

CAMRAD- 

Test 


4.34 

4.46337 

4.51864 

4.51473 

+ 4.02% 

+ 6.68% 


12.55 

12.91 

12.9917 

12.9731 

+ 3.37% 

+ 1.59% 

l 3rd 

24.99 

25.9137 

26.249 

26.2149 

+ 4.90% 

+ 2.00% 

fi 4th 

41.63 

43.3187 

43.9606 

43.8759 

+ 5.39% 

+ 1.65% 

1 5th 

63.71 

65.0184 

65.7907 

65.6459 

+ 3.04% 

- 9.35% 


It can easily be seen, that the CAMRAD-Program gives us a fairly good 
estimate of the natural frequencies that can be expected. The error varies, and 
CAMRAD fails to predict accurate natural frequencies for all higher modes. 

Even though our Finite Element model accounts for only out-of-plane flapping 
motion at this time, and disregards the coupling with the torsional and inplane 
vibrational modes, the results are in very good agreement with the measured 
data. Inclu din g these additional degrees of freedom (or allowing for coupling of 
the modes) would tend to "soften" the rotor blade and bring down the natural 
frequencies, and thus, the approach for the measured natural frequencies would 

be even better. 

Of especially great interest would be the inclusion of the torsional modes, 
since they would be coupled most with the flapping mode due to the twist and 
swept tip of the blade. In addition, the 1 st chord wise natural frequency lies close 
to the 3rd flapping natural frequency which would also tend to lead to some 

coupling. 

To judge how exact the program would calculate the natural frequencies 
given an exact value for the first natural frequency, the modulus of elasticity was 
adjusted to supply exactly the first measured natural frequency. Since (0 n is 

proportional to the square root of E, (2.3.5-1) we can simply set 

fWlcalc.'j fEJj 0,5 (2.4. 3. 1-1) 

latest J V Ead )J 
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and we receive an adjusted E as 

E a dj = 0.628383 10™ = 0.940693 10* (2.4.3.1-2) 

Using this adjusted value to simulate a blade less stiffened through modal 
coupling we find the results of Table 2.4.3.1-3. 

Table 2.4.3.1-3 Natural Frequencies Calculated with an Adjusted Modulus of 

Elasticity 


Mode 


DS3 
^mgen^ 
^en \ 


Shake Test 
[Hz] 

Change: 

DS3(FEM)- 

Test 

1st 

4.34 

4.34 


2nd 

12.47 

12.55 

- 0.59% 

3rd 

25.23 

24.99 

+ 0.989% 

4th 

42.39 

41.63 

+ 1.84% 

5th 

63.80 

63.71 

+ 0.146% 


In Figs. 2.4.3.1-3a to 2.4.3.1-3g we introduce the calculated mode shapes of 
the free-free UH-60A rotor blade along with the normalized modal residues of 
the shake test. We have used the coordinate system from the shake test used in 
the free-free analysis (i.e. going from tip to root of the blade, shown in Fig. 2.4.3. 1- 
2) to ease comparison with the test results, but in all other calculations 
(cantilevered and hinged) we have used the more common coordinate system 
going from the root of the blade outward. The measured mode shapes (residues) 
of the shake test were normalized not as usual with the maximum deflection at 
the beginning or end of the blade, but in such a manner that they coincided at 
some point along the blade with the calculated mode shapes. This was done, 
since the modal residues were assumed to vary linearly across the chord of the 
blade and calculated at quarter chord. The tip of the blade is swept back, and the 
shaker was attached to the root, and so the quarter-chord modal residues 
calculated at these points that show the maximum deflection are not considered 
to be very reliable. These points of coincidence were 4.165 m for the 1st, 3rd and 
5th mode , 5.34 m for the 2nd mode and 6.125 m for the 5th mode. 


53 




Normalized Deflection 
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UH-60A 

70 Finite Elements 



012345678 
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Fig. 2.4.3.1-3a Comparison of Calculated and Measured Flapwise Mode Shapes 


Second Mode Shape 
UH-60A 






Third Mode Shape 
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Fig. 2.4.3. 1 -3c Comparison of Calculated and Measured Flapwise Mode Shapes 
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Fig. 2.4.3.1-3d Comparison of Calculated and Measured Flapwise Mode Shapes 
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Fifth Mode Shape 
UH-60A 



Fig. 2.4.3. l-3e Comparison of Calculated and Measured Flapwise Mode Shapes 


Sixth Mode Shape 
UH-60A 
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Fig. 2.4.3.1-3f Comparison of Calculated and Measured Flapwise Mode Shapes 
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Seventh Mode Shape 
UH-60A 



Fig. 2.4.3.1-3g Comparison of Calculated and Measured Flapwise Mode Shapes 
The Calculated Data shows a high correlation with the measured mode shapes. 


The same comparison is now done for the calculated and measured natural 
frequencies for the first two inplane (chordwise) modes. Table 2.4.3. 1*4 gives the 
calculated results along with the deviations from the shake test and the 
calculated mode shapes are shown versus the measured mode shapes in Figs. 
2.4.3.1-4a and 2.4.3.1-4b. 


Table 2.4.3.1-4 Results for the First Two Inplane Modes 


Mode 

Shake Test 

FEM 

Change: 


(Hz) 

(Hz) 

FEM-Test 

1st 

25.38 

26.28 

+ 3.5 % 

2nd 

67.38 

73.27 

+ 8.7 % 
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1st Chordwise Mode 
UH-60A 










2A.3.2 Results from other Boundary Conditions 

Using DS3, the natural mode shapes and frequencies for other boundary 
conditions in a non-rotating state were calculated and are introduced in Table 
2.4.3.2-1 and Figs. 2.4.3.2-1 and 2.4.3.2-2 


Tabl^j : 4 : 3 : 2-l_jjorvRo Results for DS3 


Mode 

Hinged-Free [Hz] 

■O 

Cantilevered [Hz] 

1st 

rigid body 

2.52586 

2nd 

10.889 

17.2904 

3rd 

40.6414 

48.07 

1 4th 

75.1265 

93.7772 

l|5th 

128.8254 

152.615 


Cantileverd 
UH-60 Blade 


First four modes 



Figs. 2.4.3.2-1 First Four Modes, Non-Rotating, Cantilevered Blade 
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Hinged-Free 
UH-60 Blade 


First four modes 



Figs. 2.4.3. 2-2 First Four Modes, Non-Rotating, Hinged Blade 


2.4.4 Rotating Results 


2.4.4.1 Increase of the Natural Frequencies 

Having ascertained a high level of accuracy and convergence of our finite 
element model, we will now use this model to predict the rotating natural 
frequencies for the UH-60A Blackhawk Rotor Blade for the cases of hinged-free 

and cantilevered boundary conditions. 

The relationship between the speed of rotation and the increase of the 
natural frequencies will again be displayed in the familiar fashion using the 
dimensionless parameters 





( 2 . 33 - 7 ) 


<D= | 

Figs. 2.4.4. 

frequency 

increases. 


'grot f 

l^nij 

1-1 and 2.4.4.1-2 show the increase of the first 


(2.3.4-1) 


and second natural 


for the cantilevered and hinged rotor blade as the speed of rotation 


61 





2AA.2 Comparison between Rotating and Non-Rotating Predicted 
Mode Shapes 

Having verified that our non rotating calculated mode shapes show high 
correlation with the measured residues, we can calculate the rotating mode 
shapes and examine the change between the non rotating and rotating mode 
shapes. Figs. 2.4.4.2-1 to 2.4.4.2-4 show the change of the cantilevered natural 
mode shape of the UH-60A Rotor Blade in a non-rotating and rotating state, and 
Figs. 2.4.4.2-5 to 2.4.4.2-8 show the same for the hinged Blade. The rotational 
parameters were set to ¥ = 12 for the cantilever blade and to ¥ =1.6 for the 
hinged blade. The abscissa shows the non-normalized blade coordinate in 

meters. 

Cantilevered 


1st Mode 



Fig. 2.4.4.2-1 Rotating vs. Non-Rotating Mode Shape 
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3. Spectral Analysis of the UH 60A Flight Records 
3.1 Introduction 

After obtaining an estimate of the rotating natural frequencies and modes 
of the Blackhawk helicopter rotor blade, a first attempt was made to detect these 
natural frequencies in the flight test records. This flight test data was obtained by 
computer down-link from the Tilt Rotor Engineering Data System (TRENDS), a 
database of the NASA-Ames Research Center. The highly instrumented Flight 
Test Phase 1 (Ref. 3.1) of the Blackhawk Program provided a multitude of 
dynamic and vibrational data that was available for downloading. The data 
pertaining to the vibrational and otor characteristics available at TRENDS is 
listed in Table 3.1-2. 

The steady-state response of the rotor system will be at frequencies 
corresponding to the period of the aerodynamic driving forces. These non-linear 
excitation functions have frequencies that are integer multiples of the rotational 
speed of the rotor blade (^ , ^ ... = IP, 2 P, ...). The resulting spectral 

densities of the histories of the responses to this aerodynamic forcing function 
will be composed mainly of frequency contributions having the periodicity of IP , 
2P, ... , leaving the relatively small contributions of the transient natural modes 
"buried" in the data. During a control input of the pilot, it is hoped that these 
transient responses are excited to a measurable degree. So what we would need 
in order to uncover transients hidden in the histories is a flight test that 
contained well documented control inputs along with the response of the rotor 
blade to these input. It should be stressed at this point that this analysis does not 
constitute an in-depth evaluation of the existing database, but it was hoped that 
by a careful selection of available flights and vibrational parameters, those 
datasets would be found that showed the greatest promise of disclosing the 
desired information. The results from this first examination of the database may 
also show future researchers how to better choose and find relevant datasets in 
the database. 

Scanning the project database (See Ref 3.-2) for descriptions of the different 
flight tests performed, flight 31 was chosen for further evaluation since this flight 
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was conducted to show effects of sinusoid control inputs. The flightlog is shown 
in table 3.1-1 below. 

Table 3.1-1 : Log, Flight 31 


AIRCRAFT: 748 SINE SWEEP CONTROL INPUTS T/O GW: 18224 

FLIGHT: 31 LOCATION: AEFA CG: 3625 

FLT DATE: 22 MAY 87 COUNTERS: 3101- 3130 HRS TO INSP: 0.0 
DIRECTOR: BUCK PILOTS: CASON AND WEBRE FLT TIME: 1.2 

FLIGHT INFO: FLIGHT 3 1 WAS A SPECIAL FLIGHT THAT WAS NOT PART OF THE 
NASA PHASE I TEST PLAN. THIS FLIGHT CONSISTED OF 
SINUSOIDAL CONTROL INPUTS THROUGH ALL AXES (LONG., LAT., 
PEDAL, AND COLLECTIVE), AT TWO AIRSPEEDS: HOVER AND 120 
KIASB. 

CONFIGURATION: 

LASSIE AND BALLAST CART INSTALLED 

FLIGHT NOTES: PARAMETERS KNOWN NOT TO BE WORKING: 

ALL PARAMETERS WERE THOUGHT TO BE FUNCTIONING. 

COUNTER TYPES: MULTI AXIS SINE SWEEP CONTROL INPUTS. 

ANALOG TAPES: ITH748031.DAT 

DIGITAL TAPES: ITH748031.01X THROUGH ITH748031.30X 
VAX DATASETS: YES PLEASE. 


The available vibrational and rotor parameters are shown in Table 3.1-2 , 
along with their units. AZMRT, AZMRR, MRBR5, MRBR6 and MRBR7 were 
examined and finally AZMRT was chosen for the spectral analysis, since it 
showed the greatest change in magnitude and intensity after a control input, 
even though this parameter would be expected to show the greatest amount of 
noise. An additional factor in the selection of AZMRT was that it showed more 
intense and sharper peaks at multiples of the vibrational speed all the way up to 
10P, whereas the other vibrational parameters only showed distinct peaks up to 
5P. This is demonstrated in Fig. 3.1-1. by showing the spectral density function 
available directly from TRENDS for the vibrational parameters. 
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Table 3.1-2 Avalible Parameters 


| Rotor Parameters 


Item 

Description 

Units 

Sampling 




Freq. [Hz] 

AXMRT 

Tip accel. edgewise 

Gs 

517/1 

AZMRR 

Root acceleration flapping 

G's 

516/1 

AZMRT 

Tip accel. flapping 

G's 

516/1 

MRALSS 

MR link load aft 

Lbs 


MRBR5 

MR rear bending 50% radius 

PSI 

516/1 

MRBR6 

MRrearbending 60% radius 

PSI 

516/1 

MRBR7 

MRrearbending 70% radius 

PSI 

516/1 

MREB5 

MR edgewise bending 50% rad. 

IN-LB 

516/1 

MREB7 

MR edgewise bending 70% rad. 

IN-LB 

516/1 

MREBX1 

MR root edgewise 

IN-LB 


MRFLAP 

MR flapping 

Deg 

516/1 

MRFLSS 

MR link load forward 

Lbs 


MRLAG 

MR lead-lag 

Degs 

516/1 

MRLSS 

MR link load lateral 

Lbs 


MRNB5 

MR normal 50% radius 

IN-LB 

516/1 

MRNB6 

MR normal bend. 60% radius 

IN-LB 

516/1 

MRNB7 

MR normal bend. 70% radius 

IN-LB 

516/1 

MRNBX1 

MR root normal bending 

IN-LB 

516/1 

MRPITCH 

MR pitch 

Degs 

516/1 

MRPR 

MR pushrod load 

Lbs 

516/1 

MRSEBL 

MR shaft bending 

In-lb 

516/1 


Aircraft Parameters 


Item 

Description 

Units 

Freq 

COLLSTK 

Control position, collective 

Inches 

32/1 

LATSTK 

Control position, lateral 

Inches 

32/1 

LONGSTK 

Control position, longitudinal 

Inches 

32/1 
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FREQUENC' 


Fig 3.1-1 Spectral Density Distri 









Right 31 is subdivided into 14 datasets of about 23 second each numbered 
with counters 04-18. These counters were scanned and counters 08 and 10 were 
chosen for further evaluation. Fig 3.1-2 and Fig 3.1-3 show excerpts of these 
counters for the parameters AZMRT, the main control inputs (COLLSTK 
control position, collective and LONGSTK - Control position, longitudinal), as 
well as the rotational speed of the rotor, RPMMR. The units for these parameters 
can be found in Table 3.1-2. 
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w 


3 < 

TIME IN SECONDS 

Fig. 3.1-2 Counter 08, Flight 31 (3108) 
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3.2 Projected Natural Frequencies 


Using the Finite Element Program outlined in the previous chapters we 
calculate the first and Second rotating natural frequencies of the rotor blade at the 
given rotational speed interval to be : 


Table 3.2-1 Rotating Natural Frequencies of the Rotor Blade 


Counter 

Rot. Speed 
(RPM) 

®nl 

(Hz) 

“n2 

(Hz) 

3108 

258.0 - 261.0 

15.95 - 16.05 

40.79 - 40.89 

3110 

258.5 - 261.5 

15.97-16.07 

40.81 - 40.93 


It would of course be impossible to measure the natural frequencies as exact as 
they are given in Table 3.2-1, the main objective is to show that the natural 
frequency does not perturb much as the rotational speed varies. 

3.3 Analysis Procedure 

3.3.1 Obtaining the Spectral Density Distribution 

Power spectral densities (PSD's) of various intervals of the response histories and 

their ratios were calculated (See also Ref 3.-3, 3.-4). 

Power spectral density distributions or spectras are defined by the Fourier 
Transforms of their correlation functions: 

+ ©o 

Sxx(<0) - I RxxW e- iut <*t < 3 - 3w > 

-oo 

However, utilizing the discrete fourier transform (DFT) it is not necessary to 
calculate the explicit term for the correlation function, since the spectral density 
can be expressed directly by multiplying the fourier transform by its complex 
conjugate. This can be shown by considering the discrete formulation of the 
spectral density function, where Sfc is the DFT of the autocorrelation function 

R r . 3 


3 A discrete series can always be written as an 

~ 2irkt 27ik t, 

x(t) = a G + X a k cos(-^) + b k sin(-^r-) 

k=l 


addition of sine and cosine terms: 


7 4 



(33.1-2) 


N-l 



r=0 


R r exp 



R r is an estimate of the correlation function defined by 
N-l 

R r =rj X x s x s+r (r=0,l/2,—, N-l) 


s=o 


Inserting (33.1-3) into (33.1-1) and rearranging terms leads to 
1 « f. 2nks\ ( . 2Ttk(s-fr)^ 

Sk = 7T2 X X ^ expl 1 n J Xs+r exp l _1 N J 

1 r=o s=o x ' v ' 


(33.1-3) 


(33.1-4) 


The terms with the variable integer r can be grouped together and we can rewrite 
the expression as 



and introducing the new variable t=(s+r) 



(33.1-5) 


(33.1-6) 


If we assume x r to be periodic with period N (an assumption that has to be 
made in order to apply the DFT) , then xn+s =x s* Thus the second term of 
(33.1-6) is simply the DFT of x(s), X(s). The first term is similar to X S/ with the 
exception of the sign of the exponent. This represents X^ the complex conjugate 
of X s . This shows that the PSD can be simply calculated from the DFT by 


where 


i/ 


2jckt, 


3 k {b k }= i I x(t) cos{sin}(— j— ) dt or using complex notation for the Fourier 


Transform , X k =a k +i b k , it can be written as: 


C _ T J > 


Xk= J X(t) exp(-i(^Y^)) dt 
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S k = x s X s 


(3.3.1 -7) 


The DFT can easily be calculated by existing code for the Fast Fourier 
Transformation (FFT). 

3.3.2 Ideal Boxcar Filters in the Frequency Domain 

Ideal Boxcar filters of varying frequency bands were applied to the spectra 
obtained by the FFT, in order to supress the steady state excitation of IP, 2P and 
up. This approach did not lead to any practical results, because of the smearing of 
the peaks and the noise level inherent in the data. This analysis technique was 
found to be not applicable and quickly discarded. 

3.3.3 Mapping to the Z-Plane 

We are dealing with a discrete time series when we are evaluating the 
measurement data from the flight test. The continuous time (CT) function s(t) is 
now transformed into a discrete time (DT) function, s(nT) with values defined 
only for t =n*T , where n is an integer and T is the sampling rate. CT signals can 
be transformed to the (complex) frequency domain by the Fourier and Laplace 
transforms The result of the transformation usually provides us with additional 
insight into the operation of the system. For some systems that do not meet the 
conditions required by the Fourier transform, the Laplace transform is used, 
which involves the transformation into the s= o + to> plane. The Fourier 
transform is a special case of the Laplace transform for <7=0, that is, for input 
systems that can be represented as a superposition of sinusoidal waveforms. DT 
signals can also be represented by a Fourier transform, but instead of using the 

Laplace transform, the z-transform is used. 

We can introduce the z transform by setting up the relation between the 

the complex s plane and the complex z plane 

z = e <sT) (3.3.3-1) 

This function is a mapping of points in the s plane to points in the z plane. In 
circuit and control theory, values of s that cause a system transfer function to 
become zero or infinity (poles) provide information about the systems response 
to signals with a given frequency value. The mapping of the s plane into the z 
plane is as follows. The ico axis of the s plane (o=0) corresponds to z=exp(-io>T). 
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This is a circle of unit radius in the z plane. As ©T varies from -n to +n , this 
generates a circular path in the z plane from z=-l (/-180° ) to z— -1(/180 ). Thus the 
z transform takes a strip in the s plane between (^ ) < (0 < Of) and maps it into 

the unit circle of the z plane. The left-hand side of the s plane is mapped into the 
unit circle (s<0, 1 z I <1) . The right hand side of the s plane is mapped outside of 
the unit circle (s>0, 1 z I >1). If an analog system has poles only in the left-hand s 
plane, it is stable, and its poles will map inside the unit circle in the z plane as 
well. This mapping is shown in Fig 3.3.3-1. 



3.3.4 Digital Filter Design 

In the second part of the spectral analysis of the flight data, digital filters 
are used to filter out the high energy spectras of the first few excitation 
frequencies. The design of such a digital filter (Ref 3.-5) begins with the 
specifications of an analog filter. The parameters for a filter are the cutoff 
frequency f c , the stopband frequency f s and the characteristics of the stop and 

passbands. These characteristics are used to generate the analog transfer function, 
H(s), which is then converted to a digital filter transfer function, H(z). This 
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conversion is done by a bilinear transformation, (BLT) . Once the analog filter 
characteristics have been defined, the BLT is initiated by 


s = K 


l-Z' 1 

1+Z' 1 


(3.3.4-1) 


where s =o + tco . Depending on the manner in which the BLT is obtained , K 
can assume the values of 1, ^ or | • We will use ^ . Solving (3.3.4-1) for z leads to 

1+ctT _ = (1+oT) + itoT (3.3.4-2) 

z ” 1-oT (1-oT) - itoT 


The phase, <>, and magnitude, I z I , of z are 


izl 


(1+oT) 2 + (coT) 2 
(1-oT) 2 + (coT) 2 . 


0.5 


(3.3.4-3) 


and 


<j>(z) = tan -1 


2coT 

l-(oT) 2 -(coT) 2 


(33.4-4) 


The infinite-length i axis is mapped (nonlinearly) onto the unit circle, 
introducing a distortion of the frequency response. The distortion is in form of a 
relocation of the relevant frequencies: cutoff and stopbands. The general shape of 
the response is essentially unaffected, that is an analog lowpass filter will remain 
a digital lowpass filter, but with different cutoff and stopband frequencies. This 
shifting of the frequencies is shown in Fig. 3.3.4-1. The shifting of the frequencies 
during the transformation can be corrected by designing the analog filter to 
account for the warping due to the transformation. Then, when transformed, the 
frequencies will be at their desired locations. This prewarping can be 
accomplished by comparing the imaginary parts of the following equation that 
shows the relationship of co= 27 t f in the s plane and the angle of z, 2k Q in the z 
plane, by evaluating z on the unit circle ( z = exp(i 2 ji Q) ) and inserting into 

(3.3.4-1): 

l -[.q-UkD. 1 gtitQ _ g-ntfl 

s = T i +e -/2jtQ = T + e~» tQ 


1 tsin ttQ 

T cos jiQ 


^ i tan nQ 


= o +ko (33.4-5) 
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Comparing the imaginary parts gives us this nonlinear relation as 

(0 = y tan (ji£ 2) (3.3.4-6a) 

or solved for Q : 

£2 = — arctan (co T) (3.3.4-6b) 

7t 

Thus to design a filter suitable for filtering DT data, the frequency component co 
that you wish to filter must be expressed in the z domain variable £2. 



Fig. 3.3.4-1 Frequency Warping during the Z-Transform 


3.3.5 Power Spectral Ratios 


A section of the history of AZMRT that had little or no control inputs 
immediately before or during its interval was chosen to represent a steady state 
response. A spectral density of this steady state response was calculated and 
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related to an interval of the history that had a control input present or active, 
and the ratio of these were calculated. The sampling frequency, fsamp /°f data 

was 516.669 Hz, and so the spectra that could be examined was limited by the 
Nyquist frequency, f n = fsa ™ P ~ = 250 Hz. This limitation was of no great 

importance, since the frequency band of interest was the interval from 0 to 50 Hz, 
containing the first two natural frequencies of the rotor system that we hoped 
were excited through a control input. 

The first counter evaluated was 3110 and the time interval spanning 
0 to 3.32 seconds (1200 datapoints) was taken to be the steady-state response. The 
history is shown in Fig.3.3.5-la. Its spectral density is shown in Fig. 3.3.5-lb on a 
linear scale and in Fig. 3.3.5-lc on a logarithmic scale. The units for a power 
spectral density are [mean square/unit of frequency]. 




Frequency (Hz) 

Fig 3.3.5-la. Power Spectral Density of the Steady State Response (Ctr.3110) 
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Four different subsets ( a) to d) ) of counter 3110 were selected that showed 
promise of containing some transient response to the control inputs. The 
datasets examined are depicted in Fig. 3.3.5-2. These datasets were spectrally 
analyzed and the ratio of these to the steady-state response were calculated. Fig. 
3.3.5-3 shows the PSD of dataset a) on a linear scale, while Fig. 3.3.5-4 shows the 
same PSD on a logarithmic scale. Calculating the ratios of these datasets in a 
linear representation to the steady-state response is of course equivalent to a 
subtraction of these spectras on a logarithmic scale. 

Fig. 3.3.5-5 to 3.3.5-9 show the ratios of the response spectras of the subsets a) to 
d) of counter 3110, with the ordinate showing the magnitude of increase in 
relation to the steady state response plotted against the frequency . 
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3 20 30 40 50 6< 

Frequency (Hz) 

Fig. 33.5-6 Ratio of the PSD of b) to the Steady State Response (3110) 



Frequency (Hz) 

Fig. 33.5-7 Ratio of the PSD of c) to the Steady State Response (3110) 



Frequency (Hz) 

Fig. 33.5-8 Ratio of the PSD of d) to the Steady State Response (3110) 





Figs. 33.5-5 to 33.5-7 all showed two peaks in the ratios of the PSD in a 
frequency ranges of 17.2-17.4 Hz and 36.5-36.8 Hz that were on the order of 600 to 
5000 greater than at steady state. Increases at IP, 2P and higher were also found 
due to an increase of energy at these intervals because of an increase in the 
aerodynamic forcing functions. 17.3 Hz is an integer multiple of the rotational 
speed (4P) but the increase in the PSD is far greater than the increases at IP, 2P 
and 3P. The first natural mode was predicted at about 16 Hz, and it can be 
assumed that these peaks represent the transient responses of the excited first 
natural mode superposed with the expected excitation at 4P. The fact that the 4P 
excitation should be smaller in magnitude than the IP, 2P and 3P excitation 
supports this observation. Fig.33.5-8 is missing this characteristic peak and the 
histories responsible for these PSD's are seen to be a steady state responses at a 
new energy level with the transients damped out. Table 33.5-1 shows the 
predicted first natural frequency of the rotor blade and the frequencies at which 
the peaks in the ratio of the PSD's were located along with the deviation. 


Location of the First Natural Frequency 


Predicted Nat. 

Dataset 

Location of Peak 

Deviation 

Freq. [Hz] 


[Hz] 

[%] 

(l st /2 nd ) 


(ist/ 2 nd ) 

(1st / 2nd) 

16 / 40.9 

a) 

17.22 / 36.60 

7.60 / -10.51 

16 / 40.9 

b) 

17.23 / 35.91 

7.70 / -13.19 

16 / 40.9 

c) 

17.19 / 35.70 

7.44 / -12.71 

16 / 40.9 

d) 

17.10 / 36.60 

8.87 / -10.51 


The deviations remain at a more or less constant level which would point to a 
constant error in the rotating mode predictions. 

The next counter evaluated was number 3108. The analysis technique was 
the same as for the previous counter. In an attempt to find and identify the 
second natural mode predicted at 40.9 Hz in this counter and to better examine 
the frequency band above 30 HZ, the frequency components at IP, 2P... were 
removed by applying digital filters to the data. 

The filtering scheme involved a Butterworth lowpass filter at 3 Hz and a 
highpass at 50 Hz, the other filters were designed as stopband filters of varying 
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bandwidths and intensities of the Butterworth (BW) and Chebycheff (CH) types. 
The combination of filters that showed the best results are shown in Table 33.5-2 

and Fig 3.3.5-10. 


Table 33.5-2 Filter Combination 


Filter Type 

Order 

Bandwidth 

[Hz] 

Passband 

Ripple 

Times 

Filterd 

Gain 

[dB] 

Highp. BW 

6 

0-3 

- 

1 

- 

BW 

3 

2.04 - 4.04 

- 

2 

-16.8 

BW 

5 

6.83 - 10.83 

- 

2 

-32.9 

CH 

6 

12.12-14.12 

0.5 

3 

-17.5 

BW 

5 

16.90 - 18.40 

- 

3 

-9.3 

CH 

8 

19.82 - 23.82 

0.5 

2 

-29.8 

CH 

8 

24.11 - 28.11 

0.5 

2 

-29.1 

BW 

6 

37.1 - 41.1 

- 

1 

-9.4 

Lowp. CH 

10 

50 - oo 

0.1 

1 




This time a larger time interval from 0 to 3.87 sec (2000 datapoints), was 
chosen to represent the steady-state response, since the second natural mode 
would tend to be damped out fairly quickly. Examining longer histories increases 
the chances of having some transient response of the second mode present in the 
data. The interval shown in Fig. 33.5-1 la and Fig. 33.5-1 lb were filtered and then 
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examined for transient responses. This steady-state response is shown in Fig. 

3.3.5- 12a. Its spectral density is shown in Fig. 3.3.5-12b on a linear scale and in Fig. 

3.3.5- 12c on a logarithmic scale. 

Again, the ratio of the transient response to the steady-state was calculated 
and plotted against the frequency. These ratios are shown in Figs. 3.3.5-13 to 3.3.5- 
18. 
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Flight 31, Counter 08 
AZMRT 



TIME IN SECONDS , 

UH-60A fi/C 748 PHASE I TESTS FieC45-lla Examined Datasets 

FLT 31*8U)E SHEEP CONTROL INPUTS Wg V..M na nxammeu u sen. 

CIR 3108: myER,LHT STICK SINE: SLEEPS 




Flight 31, Counter 08 
AZMRT 



UH-60R R/C 748 PHRSE I TESTS 

FLT 31 :8INE SHEEP CONTROL INPUTS I 

CJR 3100: HOVER, LflT STICK SINE sheeps FigC.4.5-llb Examined Datasets 




Time Histo 


Time in Seconds 


Fig 3.3.5-12a. Steady State Response (Ctr.3108) 


Power spectral density 


Frequency (Hz) 


Fig 3.3.5-12b. PSD of the Filterd Steady State Response (Ctr.3108) 


Power spectral densi 



0 10 20 30 40 50 60 

Frequency (Hz) 

Fig 3.3.5-12c. PSD of the Filterd Steady State Response (Ctr.3108) 









Frequency (Hz) 

Fig. 3.3.5-13 Ratio of the PSD of a) to the Steady State Response (3108) 



Frequency (Hz) 

Fig. 3.3.5-14 Ratio of the PSD of b) to the Steady State Response (3108) 





.5-17 Ratio of the PSD of e) to the Steady State Response (3108) 
60, " T 1 


1 — js. A t U — — 

0 10 20 30 40 50 60 

Frequency (Hz) 

3 5-18 Ratio of the PSD of 0 to the Steady State Response (3108) 






Discussion of the Power Spectral Density Ratios 

All examined spectra ratios show a peak at the frequency interval of 37.9 to 
40.7 Hz, exept fot the data subset f) that had a peak at 43.9 Hz. The magnitude 
increase is less than the increase for the previous dataset, since the l ter a e 
37-41 Hz band damped both the steady state and the transient responses with a 
gain of -9.4 dB. The increase turned out to be on the order of 15 to 100 oppose to 
the increase at counter 3110 of 600-5000. Again, the second natural frequency was 
expected to be at 40.9 Hz, and it can be assumed that the peaks visible in Figs. 
3.3.5-13 to 3.3.5-18 are a result of this natural mode being excited. Ta e . . 
shows the found and expected second natural frequency. 



The results of the spectral energy analysis show great promise in being 
used to identify the transients of the natural modes buried in the response 
histories and warrant further examination. This analysis showed a definite 
presence of the first natural mode at about 17.2 Hz . Hie second mode was placed 
at 36.6 Hz for counter 3110 and at 37.9 Hz for counter 3108. Examining further 
counters in future research would lead to even better results, especially if the 
counters were cross-compared in some way. Another possibility would be to 
construct a transfer function for the assumed steady state response and one for 
the transient response and to draw conclusions from the different numerator 
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polynomials, since the denominator polynomials would have to be the same 
for the same system. 

A reason that the predictions of the rotating natural modes are too high in 
general could be explained by Fig. 3.3.6-2. As the damping increases, the natural 
modes shift slightly to the lower end of the spectrum. Our analysis neg ects 
damping effects, but in the case of the flight test measurements, aerodynamic as 
well as structural damping is present and may account for a phase sift. 


Magnitude response 



CD 

Non-dimensional natural frequency 

UJfl 

Fig. 3 .3.6-2 Magnitude of Frequency - Response as a Function of Damping 
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Appendix A 

Tables of the First Eigenvectors of Beams with Different Boundary 
Conditions and Geometric Properties 
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Cantilevered Constant Beam 



Hinged Constant Beam 


Third Mode I Fourth Mode 


0.00000 

0.0666667 

0.133333 

0.200000 

0.266667 

0.333333 

0.400000 

0.466667 

0.533333 

0.600000 

0.666667 

0.733333 

0.800000 

0.866667 

0.933333 

1.00000 


0.00000 

0.0666667 

0.133333 

0.200000 

0.266667 

0.333333 

0.400000 

0.466667 

0.533333 

0.600000 

0.666667 

0.733333 

0.800000 

0.866667 

0.933333 

1.00000 


0.00000 

- 0.176168 

- 0.342705 

- 0.484881 

- 0.590631 

- 0.651757 

- 0.663542 

- 0.624366 

- 0.535326 

- 0.399852 

- 0.223324 

- 0.0129609 

0.223915 

0.477577 

0.738882 

0.998499 


0.00000 

0.312006 

0.581464 

0.717196 

0.687227 

0.507443 

0.224895 

- 0.0968230 

- 0.390120 

- 0.593871 

- 0.661815 

- 0.568864 

- 0.315344 

0.0708468 

0.534149 

0.993338 


0.00000 

- 0.427397 

- 0.734409 

- 0.664977 

- 0.283391 

0.203671 

0.580202 

0.701849 

0.532819 

0.146648 

- 0.301164 

- 0.620699 

- 0.647359 

- 0.305840 

0.329654 

0.986899 






Cantilevered Tapered Beam 



Coor 



0.0 

0.1 

0.2 

0.3 

0.4 

0.5 

0.6 

0.7 

0.8 

0.9 

1.0 

0.00154295 

0.0147688 

0.0577830 

0.126396 

0.216554 

0.324336 

0.445950 

0.577741 

0.716185 

0.857889 

1.00000 

-0.0156330 

-0.0360000 

-0.125706 

-0.220433 

-0.272842 

-0.251821 

-0.141774 

0.0581017 

0.334526 

0.660860 

1.00000 

0.00672200 

0.0565879 

0.174609 

0.201682 

0.0921422 

-0.0906786 

-0.232354 

-0.218733 

0.0192293 

0.467178 

1.00000 


Hinged Tapered Beam 


Coord 



First Mode 
Rigid Body 


Second Mode I Third Mode I Fourth Mode" 


1 0.00000 
0.066667 
0.133333 
0.200000 
0.266667 
0.333333 
0.400000 
1 0.466667 
0.533333 
0.600000 
0.666667 
0.733333 
0.800000 
0.866667 
0.933333 
1.00000 


0.00000 

0.0666667 

0.133333 

0.200000 

0.266667 

0.333333 

0.400000 

0.466667 

0.533333 

0.600000 

0.666667 

0.733333 

0.800000 

0.866667 

0.933333 

1.00000 


0.00000 

-0.0567526 

-0.108299 

-0.144760 

-0.159416 

-0.148481 

-0.110661 

-0.0466947 

0.0410851 

0.149206 

0.273503 

0.409567 

0.553193 

0.700830 

0.850027 

1.00000 


0.00000 

0.0582599 

0.102849 

0.103761 

0.0594524 

-0.0140151 

-0.0921840 

-0.149522 

-0.164856 

-0.124921 

-0.0260174 

0.126213 

0.318904 

0.539640 

0.763966 

1.00000 


0.00000 

-0.0548090 

-0.0847653 

-0.0385510 

0.0438681 

0.102800 

0.0993520 

0.0312482 

-0.0696861 

-0.151979 

-0.164524 

-0.0757709 

0.114303 

0.375812 

0.670562 

1.00000 
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Appendix B 

Measured Residuals of the NASA - Shake Test 
(1st 5 Flapping Modes) 
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1st Fla 


lit 


1.046 

870.050m 

864.000m 

340.000m 

259.000m 

-113.000m 

-259.110m 

-468.700m 

-587.690m 

-932.680m 

-986.740m 

-1370 

-1.164 

-1336 

-1.407 

-1.409 

-1.361 


2nd Fla 


-610.000m 

-1.413 

440500m 

-262500m 

1.339 

661.246m 

2.256 

1.392 

2.500 

1.857 

2.490 

1.718 

2.210 

1.485 

1.420 

766500m 

600.000m 

-115.430m 


>Tr 




Acccloro- 
|| meter 
position 

Residue 

(m=milli) 



Acceloro- 

mcter 

position 

1 

4.040 

21 1 

2 

3.530 

22 3 

3 

976.000 m 

23 : 

4 

726.000 m 

24 

5 

-710.700 m 

25 

6 

-802.000 m 

26 

7 

-1.864 

27 

8 

-1.940 

28 

9 

-2340 

29 

10 | 

-2.153 ! 

30 

11 

-1.736 

31 

12 

-1.605 

32 

13 

-705.000 m 

33 

14 

-407.000 m 

34 

15 

614.000 m 

35 

16 

926.000 m 

36 

17 

1.589 

37 

18 

1.982 

38 

19 

2.175 


20 

2.525 



Residue 

(m=milli) 


-3.217 
-1.772 
-3389 
-1.710 
-2345 
-913.000 m 
-895.000 m 
572.000 m 
741.7000 m 
-2.132 


Acccloro- 

mctcr 

position 

Residue 

(m=milli) 

Acceloro- 

mcter 

position 

Residue 

(m=milli) 

1 

3.140 

21 | 

-2306 

2 

2.356 

22 

-700.000m 

3 

1.660 

23 

-2.219 

4 

527.000m 

24 

-936.000m 

5 

393.000m 

25 

-2.763 

6 

-1319 

26 

-420.000m 

7 

5.000m 

27 

-1362 

8 

-1.425 

28 

534.000m 

9 

416.000m 

29 

-1.013 

10 

724.000m 

30 

1.521 

11 

1.530 

31 

-706.000m 

12 

306.000m 

32 

1.960 

13 

1.602 

33 

-1 .250 

14 

1.132 

34 

1.642 

15 

1.408 

35 

-2.210 

16 

1.351 

36 

698.000m 

17 

425.000m 

37 

414.000m 

18 

874.000m 

38 


19 

-966.000m 



20 

35.500m 
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ll 5th Hap |1 

Acceloro- 

Residue 

Acceloro- 

Residue || 

meter 

position 

(m=milli) 

meter 

position 

(m=milli) 

1 

-2500 

21 

1.243 

2 

-2.983 

22 

783.000m 

3 

1.250 

23 

-1.081 

4 

-78.000m 

24 

-636.000m 

5 

3.040 

25 

-2.641 

6 

1.042 

26 

-1.497 

7 

2.880 

27 

-2510 

8 

681.000m 

28 

-1.074 

9 

1.190 

29 

-1.060 

10 

-718.000m 

30 

370.000m 

11 ! 

-408.000m 1 

31 

419.000m 

12 

-1.870 

32 

1.691 

13 

-498.000m 

33 

258.000m 

14 

-1.774 

34 

1.829 

15 

830.000m 

35 

-1.077 

16 

-559.000m 

36 

659.000m 

17 

2.349 

37 

709.000m 

18 

888.000m 

38 


19 

2.681 



20 

1.466 
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Appendix C 

Finite Element Program for Rotating Structures 
Example for Hinged UH- 60A Rotor Blade 



The Finite Element Program was written on a NeXT computer using the 
MATHEMATICA language (see Ref. 06). The authors native 
German, so some of the internal variables are 

German These are listed below to enhance understanding of the code.Th 
variables used in the code are consistent to those in the report, Greeksymbok are 
spelled out. All other variables are self-explanatory or are defined in the 
comments. Three different fonts were used in the listing to help discern between 

input, output and comments: 

Code Text 

Comment 


Output Text 


Variables: 

anzahl 

PSI 

elementlaenge 

elementlaengez 

ort 

mquerz, mquer 

ne, n 
kele, kstr 

mele, mstr 

schleife 

aktele 

zwisch, zw 

letztespalte 

kuudach 

mue 

teiler 

bild 

nust 

mnust 

lges 


parameter governing the amount of plots 

rotational parameter Y 

length of the finite element 

as above, in reference frame "z" 

location along the blade or beam 

kg 

mass distribution [ -£] 
number of elements 

part of the element (ele) and the structure (str) stiffness 
matrix 

see above, with mass matrix 
loop variable 
current element 
an intermediate result 
last column 

K uu , reduced stiffness matrix 
It 

divisor 

picture 

numerical derivative of u, the eigenvector 
matrix of the above 
total length 
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Structure of Example Program 


Preprocessor 


Input data of structure 

- divide structure 

- define element length 

- siffness and mass distribution 

- rearrange coordinates to run 
from root to tip 


^ Assembly of the System Matrices 

- define element matrices 

- assemble structure 

^ - static cond ensation 


V. 


Calculate Eigenvalues and Eigenvectors 

- Eigenvalues 

- for well conditioned matrices 

- with normalized matrices 

- Eigenvectors 

- calculation . 

- normalize and calculate the system matrix 

- plot eigenvectors 


i 


The Rotator 


A 


Rotating the Structure 

- calculate the derivatives of the eigenvectors 
. calculate the dynamic stiffness matrix 




r 


i 


Output Section 




Outpiit 

- calculate the eigenvalues of the rotating 
structure 

- find the rotating eigenvectors 

- transform back to physical coordinate 
system 

- add dynamic stiffness 

- find new system matrix 

- plot non-rotating and rotating eigenvectors 

- find polynomial approach to eigenvectors 

- plot polynomials 

- plot comparison of rotating and 
non-rotating eigenvectors 

Save Data to Optical Disc 




J 


1 04 






PREPROCESSOR BEGINNING 

Th. Breitfeld 1990 

DATA 

Parameter 

anzahl=5; 

number of modes to be considerd : modes 

modes=5 

PSI=1.6; 

Structuredata 

The Blade is divided into 3 Parts 

Part I : 0....2.0 m (nl Elements) 

Part 2 : 2 JO....6J0 m (n2 Elements) 

Part 3 : 6J0 .... 7,29 m (n3 Elements) 

Tip of Blade: xi=0 

it of Elements : n 

11=1.29 

12=4 

13=2 

Print[ll +12+13] 
nl=25 
n2=20 
n3=25 

n=nl+n2+n3 

ne=n; 

Lengths of the Single Elements 

Part 1 (tip...) 
ii=0; 

Do[ 

ii=ii+l; 

elementlaenge[ii]=ll / nl; 

/{nl}] 

Part 2 
Do[ 

ii=ii+l; 

elementlaenge[ii]=12/ n2; 

/{n2}] 

Part 3 (...root) 

Dot 

ii=ii+l; 

elementlaenge[ii]=13/ n3; 

/{n3}] 

Coordinate Table xi (Tip....Root) 

Table[xi[ii]=0 / {ii / l,n]]; 

xi[l]=Nlelementlaenge[l]]; 

ort=l; 

Dot 
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ort=ort+l; 



xi[ort] =(xi[ort-l ] +N[elemen tlaenge[ort]]) ; 


,{n-l)] 

Stiffness Tables [cm A 4] 
n=l; 

"TIP VALUE"* 

While[xi[ii] <= .1016 , {i[ii]=308 10 A (-8) , ii++}] 
While[xi[ii] <= .2972 , {i[ii]=370 10 A (-8) , ii++}] 
While[xi[ii] <= .4928 , {i[ii]=614 10 A (-8) , ii++}] 
While[xi[ii] <= .8222 , {i[ii]=826 10 A (-8) , ii++}] 
While[xi[ii] <= 1.1516 , {i[ii]=929 10 A (-8) , ii++}] 
While[xi[ii] <= 1.5555 , {i[ii]=945 10 A (-8) , ii++}] 
While[xi[ii] <= 1.9594 , [i[ii]=978 10 A (-8) , ii++}] 
While[xi[ii] <= 2.3633 , (i[ii]=978 10 A (-8) , ii++}] 
While[xi[ii] <= 2.8075 , {i[ii]=978 10 A (-8) , ii++}] 
While[xi[ii] <= 3.2517 , (i[ii]=964 10 A (-8) , ii++}] 
While[xi[ii] <= 3.6959 , (i[ii]=961 10 A (-8) , ii++}] 
While[xi[ii] <= 4.1401 , {i[ii]=961 10 A (-8) , ii++}] 
While[xi[ii] <= 4.5668 , {i[ii]=933 10 A (-8) , ii++}] 
While[xilii] <= 4.9935 , {i[ii]=924 10 A (-8) , ii++}] 
While[xi[ii] <= 5.4202 , {i[ii]=924 10 A (-8) , ii++}] 
While[xi[ii] <= 5.8469 , {i[u]=924 10 A (-8) , ii++}] 
While[xi[ii] <= 6.2736 , {i[ii]=924 10 A (-8) , ii++}] 
While[xi[ii] <= 6.9086 , {i[ii]=1149 10 A (-8) , ii++}] 
While[xi[ii] <= 7.3531 , {i[ii]=3900 10 A (-8) , ii++}] 
* ROOT VALUE”; 

Mass Tables [kg/m] 
ii=l; 

“TIP VALUE"; 

While[xi[ii] <= .1016 , (mquer[ii]=3.99 , ii++}] 
While[xi[ii] <= .2972 , [mquer[ii] =12.55 , ii++}] 
While[xi[ii] <= .4928 , (mquer[ii] =16.75, ii++}] 
While[xi[ii] <= .8222 , {mquer[ii]=17.79, ii++}] 
While[xi[ii] <= 1.1516 , [mquer[ii]=16.27, ii++}] 
While[xi[ii] <= 1.5555 , {mquer[ii]=13.03 , ii++}] 
While[xi[ii] <= 1.9594 , {mquer[ii]=13.1 , ii++}] 
While[xi[ii] <= 2.3633 , (mquer[ii]=13.17 , ii++}] 
While[xi[ii] <= 2.8075 , |mquer[ii] =11.65, ii++}] 
While[xi[ii] <= 3.2517 , {mquer[ii]=10.55, ii++}] 
While[xi[ii] <= 3.6959 , [mquer[ii] =10.41, ii++)] 
While[xi[ii] <= 4.1401 , {mquer[ii] =10.62, ii++}] 
While[xi[ii] <= 4.5668 , {mquer[ii] =10.34, ii++}] 
While[xi[ii] <= 4.9935 , {mquer[ii]=10.34, ii++}] 
While[xi[ii] <= 5.4202 , [mquer[ii] =10.34, ii++)] 
While[xi[ii] <= 5.8469 , {mquer[iij =10.20, ii++}] 
While[xi[iij <= 6.2736 , {mquer[iij=10.27, ii++}] 
While[xi[ii] <= 6.9086 , {mquer[ii]=10.34 , ii++)] 
While[xi[ii] <= 7.3531 , (mquer[ii]=30.06 , ii++}] 
"ROOT VALUE"; 


Testdata 

ro=7800 

ii=0; 

Do[ 

ii=ii+l; 

e[ii]=0.68 10 A 11; 

/{n}] 

The Coordinate Shuffler 

3d* - tip. root changes to xi - root...tip 

x -root...tip 

ii=. 

TablelxIii]=0 / {ii,l / n}]; 

ii=0; 

Do[ 

ii=ii+l; 

x[ii]=(ll+12+13)-xi[ne-ii]; 

,{n-l}] 

x[ne]= 11+12+13; 
ii=0; 


Do[ 

ii=ii+l; 

xi[ii]=x[ii]; 

Renumber the values for stiffness and mass and elementlength... 

ii=0; 

Do[ 

ii=ii+l; 

iz[ii]=i[ne+l-ii]; 

mquerz[ii] =mquer [ne+ 1 -ii] ; 

elementlaengez[ii]=elementlaenge[ne+l-ii]; 

An)] 

Back to the org. names 

ii=0; 

Do[ 

ii=ii+l; 

i[ii]=iz[ii]; 

mquer[ii] =mquerz[ii]; 
elementlaenge[ii] =elementlaengez[ii] ; 

,tn)] 

Hinged-Free 
dof=n 2 +1 


ne=n; 
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PREPROCESSOR END 



Assembly of the FEM-Stiffnessmatrix 
Element stiffness matrix 

ne= number of the element 
“ 4x4 matrix 


kel=Array[kele,{4 / 4)]; 


kele[l,l]=12 a 

kele[l/2l=6 elementlaenge[aktele]a 
kele[l,3]=-12 a 

kele[l,4]=6 elementlaenge[aktele]a 

kele[2,2]=4 elementlaenge[aktele] elementlaenge[aktele] a 
kele[2,3]=-6 elementlaenge[aktele]a 

kele[2,4]=2 elementlaenge[aktele] elementlaenge[akteleja 
kelel3 / 3]=12 a 

kele[3 / 4]=-6 elementlaenge[aktele]a 

kele[4,4]=4 elementlaenge[aktele] elementlaenge[aktele]a; 
run=l; 

Do[ 

run=run+l; 


ll=run-l; 

schleife=0; 

Do[ 

schleife=schleife+l; 

kel[[run,schleife]]=kel[[schleife / run]]; 

,{3}]; 

Assembly of the Structure 
Set kstr to zero 

kstr=Array[kst,{dof,dof}]; 

run=0; 

schleife=0; 

Do[ 

schleife=0; 

run=run+l; 

Do[ 

schleife=schleife+ 1 ; 
kstr[[run,schleife]]=0; 

,{dof}] 


,(dof}]; 

first element (hinged boundary conditions) 
flktGlG = l * 

a=elaktele] i[aktele] /(elementlaengelaktele] A 3); 


kstr[[l,l]]=kel[[2,2]] 

kstr[[l / 2]]=kel[[2 / 3]] 

kstr{[l / 3]]=kel[[2 / 4]] 

kstr[[2 / 2]]=kel[[3 / 3]] 

kstr[[2 / 3]]=kel[[3,4]] 

kstr[[3 / 3]]=kel[[4 / 4]]; 
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second until n-th element 
run=l; 


faktor=-2; 

Do[ 

run=run+l; 
aktele=run; 
faktor=faktor+l; 
index=aktele+faktor+l; 
a=e[aktele] i[aktele]/(elementlaenge[aktele] A 3); 

kstr[(index,index]]=kstr[[index / index]]+kel[[l,lj]; 

kstr[[index / index+l]]=kstr[[index4ndex+l]]+kellll,2 ; 

kstr[[index4ndex+2]]=kstr[[index,index+2]]+kel[[l / 3JJ; 

kstrl[mdex4ndex+3]]=kstr[[index4ndex+33]+kellll / 4JJ; 

kstr[[index+l / mdex+l]]=kstr[[index+l,index+l]]+kel[[2 / 2]l; 

kstr[[index+l / index+2]]=kstr[[index+l / mdex+2]]+kel[2 / 3J; 
kstr[Undex+l / index+3]]=kstr[[index+14ndex+3] +ke 2,4 ; 
kstr[[index+24ndex+2]]=kstr[[index+24ndex+2 +ke 33 J; 

kstr[[index+24ndex+3]]=kstr[[index+2,index+3]]+kel[[3,4]]; 

kstr[[index+3 / index+3lj=kstr[[index+34ndex+3J]+kelL[4 / 4JJ; 


Symmetry 

matrix=kstr; 

zwisch=Transpose[matrix]; 

matrix=matrix+zwisch; 

ii=0; 

Dot 

ii=ii+l; 

matrix[[ii4i]]=nxatrix[[ii4i]]-zwisch[[ii4ij]; 


,{dof}]; 

kstr=matrix; 

Assembly of the mass matrix 
Set mel to zero 


run=0; 

schleife=0; 

Dot 

schleife=0; 

run=run+l; 

Dot 

schleife=schleife+l; 
mel [[rurt/Schleife] ] =0; 


/{4)1 

,{ 4 ) 1 ; 

aktele~v 

mel[[l ,1 ] 1 =mquer [aktele] elementlaenge[aktele] / 2; 
melt [331] =mquer[ak tele] elementlaenge[aktele]/2; 
set mstr to zero 

mstr=Array[mst,(dof,dof)]; 
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run=0; 

schleife=0; 

Dot 

schleife=0; 

run=run+l; 

Do[ 

schleif e=schleif e+ 1 ; 
mstr[ [run,schleife]] =0; 

Adom 

,{dof)]; 

first element (Cantilever boundary conditions) 
slctdc — | * 

mstr[[l / l]]=mel[[2 / 2]] 

mstr[[2,2]]=mel[[33]]; 

second until n-th element 

run=l; 

faktor=-2; 

Do[ 

run=run+l; 

aktele=run; 

faktor=faktor+l; 

index=aktele+faktor+l; 

mstr[[index / index]]=mstr[[index4ndex]]+mel[[l / l]]; 

mstr[[index+l / index+l]]=mstr[[index+l / index+l]]+mel[[2 / 2]]; 

mstr[[index+2 / index+2]]=mstr[[index+2,index+2]]+mel[[3 / 3]]; 

mstr[[index+3 / index+3]i=mstr[[index+3 / index+3]]+mel[[4,4]]; 

,{ne-l)]; 

Static Condensation 
Massmatrix 
Rearrange the DOFs 

(Rearrange Elements Naturally And Transform Equations 

to Reduced Stiffness matrices) 

m=Array[mm / {dof,dof}]; 

run=0; 

schleife=0; 

Do[ 

schleife=0; 

run=run+l; 

Do[ 

schleife=schleife+l ; 
m[[run,schleife] ] =0; 

Adom 

A dof)]; 
colums 
run=l; 
index=l; 

Dot 

m[[index]]=mstr[[run+l]]; 

index=index+l; 

run=run+2 


1 1 0 



,{(dof+l)/2-l}]; 

rows 

run=0; 

Do[ 

run=run+l; 

m[[run,run]]=m[[run,(2 run)]]; 
m[[run,(2 run)]]=0; 

,{((dof+l)/2)-l}] 

Structurematrix 
rearrange the DOFs 
k=Array[kk,{dof / dof)]; 
run=0; 
schleife=0; 

Do[ 

schleife=0; 

run=run+l; 

Do[ 

schleife=schleife+l; 

k[[run,schleife]]=0; 

,{dof]] 

/{dof}]; 
rows, colls 

letztespal te=ks tr[ [dof] ]; 
run=0; 

Do[ 

run=run+l; 

k[[((dof+l)/2-l)+run]]=kstr[[2 run -1]]; 
k[[run]]=kstr[[(2 run)]]; 

,((dof-l)/2]]; 
k[[dof]]=letztespalte; 
kstr=Transpose[k]; 
letztespal te=kstr[ [dof] ]; 
run=0; 

Do[ 

run=run+l; 

k[[((dof+l)/2-l)+run]]=kstr[[2 run -1]]; 
k[[run]]=kstr[[(2 run)]]; 

,{(dof-l)/2}]; 
k[[dof] ] =letztespal te; 
k=Transpose[k]; 

Condensation 
Massmatrix 
remove zeros 
ii=0; 

Do[ 

ii=ii+l; 

list=m[[ii]]; 

zw=Partition[list / (dof-l)/2]; 

m[[ii]]=zw[[l]]; 
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,{(dof-l)/2}J; 

ii=.; 

m=Table[m[[ii]] / {ii / l / (dof-l)/2)]; 

Stiffness matrix 
zw=Array[zz,{dof,dof}]; 

the following variables correspond to those in chapter 2 . 1.5 in their notation 

kuu=Array[uu / {(dof-l)/2 / (dof-l)/2}]; 

kuf=Array{uf,{(dof+l)/2,(dof-l)/2}]; 

kfu=Array[fu,{(dof-l)/2,(dof+l)/2}]; 

kff=Array[ff / {(dof+l)/2,(dof+l)/2}]; 

ii=0; 

Do[ 

ii=ii+l; 

zw[[ii]]=Partition[k[[ii]] / (dof-l)/2]; 

,{dof}]; 

ii=0; 

Do[ 

ii=ii+l; 

kuu[[ii]]=zw[[ii]][[l]] ; 

/ {(dof-l)/2}]; 


ii=0; 

Do[ 


ii=ii+l; 

kuf[[ii]]=zw[[ii+(dof-l)/2]][[l]]; 

/ {(dof+l)/2}]; 

ii=0; 

Do[ 

ii=ii+l; 

rot[ii]=RotateLeft[k[[ii]],(dof-l)/2]; 
zw[ [ii]] =Partition[rot[ii],(dof+l) /2]; 

,{dof}] 

ii=0; 

Do[ 


/ {(dof-l)/2}] 

ii=0; 

Dot 


ii=ii+l; 

kfu[[u]]=zw[[u]][[l]]; 


ii=ii+l; 

kff [l a H= zw t[((dof-l)/2)+ii])[[l]]; 

/ {(dof+l)/2}] 

The reduced Stiffness matrix; 

kuudach= kuu- Transpose[kuf| . Inversefkff] . kuf; 

k=kuudach; 

dof=(dof-l)/2; 
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Search for the Eigenvalues 

Normal Apprch (well conditioned matrix) 

mm=Inverse[m]; 

eig=Eigenvalues[k . mm]; 

eig A .5; 

With the normalized matrices 

the letter "n" or the preposition "norm" depicts a normalized entity... 

matrix=. 

matrix=k 

dmatrix= matrix Transpose[matrix] 

lamdas=Eigenvalues[dmatrix] 

qnorm=Max[lamdas] 

knorm=qnorm A .5 

nk=matrix 1/knorm 

Min[Abs[Eigenvalues[nk]]] 

if this value is « 1 ,then the matrix is indeed singular 

matrix=. 

matrix=m 

dmatrix= matrix Transpose[matrix] 

lamdas=Eigenvalues[dmatrix] 

qnorm=Max[lamdas] 

mnorm=qnorm A .5 

nm=matrix 1/mnorm; 

nmm=Inverse[nm]; 

nk . nmm; 

mue=Eigen values [%]; 
lamq=mue knorm/mnorm; 

% A .5 

lower eigenvalues (roots) belong to the lower eigenfunctions 
roots=lamq; 

Reinsertion of the Eigenvalues into the matrices 
lamq=.; 

eigprob=k - lamq IdentityMatrix[dof] m; 
ii=0; 

Do[ 

ii=ii+l; 

lamq=roots[[dof+l-ii]]; 

{d,p)=Eigensystem [eigprob]; 
y=Table[0,{dof}]; 


“the lowest term of “d“ must be found, sothe correct 

position of y is set to 1 

pos=0; 

jj=i; 

tst=Abs[d]; 

Do[ 

If [ ts t [ [ j j] ] = =M i n [ ts t] ,pos = j j ,pos =pos]; 

jj=jj+h 

,{dof}]; 
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yl[pos]]=l; 
u[ii]=y . p; 

,{dof}] 

"Normalizing the eigenvectors 

run=0 

Do[ 

IfIAbs[Max[u[run]]]>Abs[Min[u[run]]] / teiler=Max[u[run]] / teiler=Min[u[run]] 

l 

u[run]=u[run] 1/teiler; 
precede the eigenvectors with a zero (for better plotting): 

ii=0; 


ii=ii+l; 

akt=u[ii]; 

uplot[ii]={0,akt); 

uplot[ii]=Flatten[uplot[ii]]; 

,{dof}]; 

Variables to this point; 

u(i) - normalised eigenvektors . . . ... v 

uplot(i) - normalised eigenvektors w/0 at beginning (for plotting) 

polyu(i) - polynominal apprch for the above 

m - kondensed mass matrix 

k -kondensed stiffness matrix 

TOO ts - squares of the eigenfrequencies (==lamq) 

dof -degrees of freedom 

n -number of finite elements(==dof) 

anzahl -plotting parameter 

Generalised Stiffness and Mass Matrix 


evek=Table[u[ii] / {ii,l / dof}]; 

U=MatrixForm[evek]; 

UNORM=Transpose[evek]; 
mgen=Transpose[UN ORM] . m . UNORM; 
kgen=Transpose[UNORM] . k . UNORM; 

omegaqu=Table[kgen([nn / nn]]/mgen[[nn,nn]] / {nn / l / n}J, 

% A (.5) 

draw the eigenvectors 
run=0; 

Dot 

run=run+l; 

"coordinate of the current element"; 

xi[l ] =elementlaenge[l ]; 

ort=l; 

Dot 

ort=ort+l; 

xi[ort] =xi[or t-1 ] +elementlaenge[ort]; 

,{n-l}]; 

datapairs=Table[coords,{n)]; 


114 



Do 


bild[run] 

,{anzahl}] 

14 
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0.64 


0.44 
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The ROTATOR 

The Derivatives of the Eigenvectors 
Procerdure for obtaining the numerical derivatives 

ii=0; 

Do[ 

ii=ii+l; 

mnust=Array[nust,{dof,dof}]; 

deltax=l / n; 

"first"; 

nust[ii,l]=u[ii][[l]] /( xi[l]/lges); 

nust[ii,dof]=(u[ii][[dof]]-u[ii][[dof-l]]) / ( xi[dof]/lges-xi[dof-l]/lges); 

"others "; 

run=l; 

Do[ 

run=run+l; ,, . x 

nust[ii,run]=(u[ii][[ r un+l]]-u[ii][[run-l]] )/ ( xi[run+l]/lges-xi[run-l]/lges); 

,{dof-2}]; 

ax[ii]=Table[nust[ii / kk] / (kk / l/dof) ]; 

,{modes}]; 

mnust=Array[ax / {modes)]; 


set gamma (nxn) -0 

schleife=0; 

run=0; 

Do[ 

schleife=schleife+l; 

run=0; 

Dot 



run=run+l; 

gamma[schleife / run]=0; 


/{n}] 


deltax=. 
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deltax[l]=xi[l]/lges; 
deltax[n]=xi[n] /lges-xi[n-l] /lges; 
ii=l; 

Do[ 

ii=ii+l; 

deltax[ii]=xi[ii+l]/lges-xi[ii] /lges; 

>-2)1 

check if sum of deltax=l 

ii=. 

SumldeltaxliiLlii/l/n}] 

0.996104252400547 

jj=* 

mm=.; 

mm=0; 

rr=0; 


Dot 


rr=rr+l; 

Print[rr]; 


mm=0; 


Do[ 


mm=mm+l; 

gamma[rr,mm] 


=Sum[m[[i41] xi[i]/lges Sum[mnust[[rr,jj]] 

mnust[[mm,jj]] deltaxljj], {jj,l,i}]/ {i,l,dof}] 


,{modes}]; 

/ {modes}]; 


1 

2 

3 

4 

5 

Grm=Array[gamma,(dof,dof}]; 


omega=.; 

kdyn=omega A 2 Grm; 


Output section 

omega=0; 
kges=kgen+kdyn; 
kges . Inversetmgen]; 

Eigenvalues[%]; 

freq=% A .5; 

freq=Sort[%] 

wnl=freq[[2]] 


68.41862788567951 


run=0; 

Do[ 

omega = (run 4/10) A .5 wnl; 

Print[N[run 4/10]]; 

kges=kgen+kdyn ; 

z=kges . Inver se[mgen]; 

tt=Eigenvalues[z]; 

xx=tt A .5; 

xx=Sort[xx]; 

Print[xx[[2]]]; 

xxx=(xx[[2]]/wnl) A 2; 

Print[xxx]; 

Print[""]; 

run=run+l; 

/{ 6 ]] 

0 . 

68.4186 

1 . 


0.4 

135.412 

3.91711 

0.8 

178.192 

6.7831 

1.2 

212.264 

9.6251 

1.6 

241.444 

12.4533 

2 . 

267.382 

15.2727 

omega=0; 
kges=kgen+kdyn; 
z=kges . Inverse[mgen]; 
Eigenvalues [z]; 


freq=% A .5; 

freq=Sort[%] 

wn2=freq[[3]] 


225.3581507539162 



run=0; 

Do[ 

omega=(run 4/10) A .5 wnl; 

Print[N[run 4/10]]; 

kges=kgen+kdyn ; 

z=kges . Inverse[mgen]; 

tt=Eigenvalues[z]; 

xx=tt A .5; 

xx=Sort[xx]; 

Print[xx[[3]]]; 

xxx=(xx[[3]] / wn2) A 2; 
Print[xxx]; 

Print[""]; 

run=run+l; 

,{ 6 }] 

0 . 

225.358 

1 . 


0.4 

297.736 

1.74549 

0.8 

354.318 

2.47195 

1.2 

402.216 

3.18545 


1.6 

444.465 

3.88981 


2 . 

482.682 

4.58749 


Find the Rotating Eigenforms 
wnl 

68.41862788567951 


Rotational Speed 
omega=PSI A .5 wnl 
86.5434794009031 
ii=. 

UNORM^ransposetTabletuIiiLUi/Udof)]]; 
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tbl =Inverse[Transpose[UNORM]l . kgen . taversetUNORM]; 
tb 2 =Inverse[Transpose[UNORM]] . kdyn . InversetUNORM], 


tbmass=Inverse[Transpose[UNORM]] . mgen . InverselUNORMl; 

c=Inverse[tbges]; 

d=c.tbmass; 

{eval,evek}=Eigensystem[d]; 


Urot=Transpose[evek]; 

”Urot=strukturmatrix (not Normalized) 
m urot[il=Eigenvektorens (rotating) 

run=0 

Do[ 

Tun=run+1; 

urot[run]=evek[[run]]; 

MatrixForm[urot[run]]; 

,{n}] 

"Normalizing" 

run=0 

Dol 

IflAbslMax[urotlrunl]]>AbstMinlurot[run]]] / teiler=Max[urot[run]],teiler=Min[urot[ 

run]]]; 

urot!run]=urot[run] / teller 

,{n}] 


"formulate the rotating structure matrix" 


run=0; 

Do[ 

run=run+l; 

zwisch=T able[urot[nn] , (nn,l ,n] ]; 


UNORM=MatrixForm[Transpose[zwisch]]; 

mgenrot=Transpose[Urot] . tbmass . Urot 

mgenrot=TableUf[mgenrot[[kk / j]]<10 A (-10) / 0 / mgenrot[[kk,j]]],(kk,n} / {j,n}] 

MatrixFormlmgenrot]; 
kgenrot=Transpose[Urot] . tbges . Urot 

kgenrot=Table[If[kgenrot[[kk / j]]<10 A (-6) / 0,kgenrot[[kk / j]]] / {kk / n} / {j, }] 


MatrixFormlkgenrot]; 


omegaqurot=Table[kgenrot[[nn,nn]] / mgenrotKn^nnllAnnJ.n}]; 

Sort[% A (-5)] 


(l/eval) A .5 



Plot nonrotating eigenvectors 
run=0; 

Do[ 

run=run+l; 

" coordinate of the current element ; 

xi[l]=elementlaenge[l]; 

ort=l; 

Dot 

ort=ort+l; 

xi[ort]=xi[ort-l]+elementlaengelortJ; 


datapairs=Table[coords / ln) J; 


ii=0; 

Dot 


datapairs[[U]]=()d[u] / utrun][[ii]]); 


bild[run]=ListPlot[datapairs]; 

dptrun]=datapairs; 
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rotating eigenvectors 

run=0; 

Do[ 

run=run+l; 

"coordinate of the current element 

xi[l]=elementlaenge[l]; 

ort=l; 

Do[ 

ort=ort+l; 

xi[ort] =xi[ort-l]+elementlaenge[ort]; 

,{n-l}]; 

datapairs=Table[coords,{n}]; 

ii=0; 

Do[ 

ii=ii+l; 

datapairs[[ii] ] = {xi [ii],urot[run] [[ii] ]}; 

bildrot[run]=ListPlot[datapairs]; 
dprot[run] =datapairs; 

,{anzahl}] 
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"polynominal approach non-rotating” 


run=0 

ii=. 


x=. 

Do[ 

run=run+l; 






" desreee of approach: Nr. of eigenform+3 

».**** ********* ; 

polyu[run]=Fit[dp[run] / Table[x A ii / {ii,0 / run+3}] / x]; 

,{anzahl}] 

run=0 


Do[ 


run=run+l; 

fkt[run] =Plot[polyu[run],{x / 0,lges} ]; 


,(anzahl}] 
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" rotating" 

run=0 

ii=. 

x=. 

Do[ 

run=run+l; 

«****»**********+*** *********************"; 
»*++* 4 *+**+++++***+++++++**++***+*+*+***+ ” ; 

polyurot[run]=Fit[dprot[run],Table[x A ii / (ii,0 / run+2}] / x]; 

,{anzahl}] 

lges 

7.29 

run=0 

Do[ 

run=run+l; 

fktrot[run]=Plot[polyurot[run],{x / 0 / lges}]; 

,{anzahl}] 
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Save data 

Write["nrothil ",u[l] ] 

Write["nrothi2",u[2]] 

Write["nrothi3",u[3]] 

Write["nrothi4" / u[4] ] 

Writer , nrothi5" / u[5]] 

Write["rothil" / urot[l]l 

Write["rothi2" / urot[2]] 

Write["rothi3" / urot[3]] 

Write["rothi4" / urot[4]] 

Write["rothi5" / urot[5]] 

ii=- 

xi70=NlTable[xi[ii] / (ii / l,n}]] 
Write["xi hin (70)",xi70] 
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Abstract 

Results from a program to instrument the 
rotor blades of a light autogiro are 
The work was initiated to provide additional 
data on rotor dynamic response as well as 
investigate practical implementation issues on 
Zu*o( blade-mounted instrumentation for 
rotor state feedback. A description of the 
aircraft and rotor electronics hardware desig 
and installation is given, along with results to 
date from the initial flight test program for 
complete system check-out. 


Introduction 

Recent attempts to expand both the 
fidelity of engineering models for r0 ^ cra 
systems and the frequency bandwidth for 
helicopter flight controllers have required 
more accurate models of coupled rotor and 
airframe dynamics. The complexity of the 
engineering modeling problem c »upled w, th 
the general lack of sufficiently detailed flight 
test data, have made improvements in 
complete rotorcraft aeroelastic pred'Ctions 
difficult. Efforts to expand the available data 
for correlation exercises are underway UL but 

such programs may experience funding and 

operational delays along with regulatory 
hurdles that preclude rapid turnaround for 
timely engineering research efforts. 

As a means of addressing this data 
deficiency, a program is underway at Pnnceton 
University to instrument a Bensen B8-M 
Gvroglider for towed-flight investigations 
(Figure 1) 15]. This aircraft is an extremely 
simple teetering rotor autogiro whose power is 
supplied by towing the vehicle behind an auto 
along a runway or suitably paved surface. The 
University's inactive runway at the Forrestal 
Campus has served as the vehicle's testing 
ground. The use of a simple test aircraft 
provides the ability to perform fundamental 

Presented at the 16th European Rotorcraft 
Forum, Glasgow, Scotland, September 18- , 

1990. 


aeroelasticity experiments on a full size 
vehicle without the additional burden of 
maintenance manpower associated with a 
production helicopter. The added capability 
of an "in-house" test vehicle affords the 
researcher the luxury of planning and executing 
tests that are driven by the nature of the test 
data and not the predetermined schedule of 
the test program. 

The desire to conduct instrumented rotor 
experiments on full size aircraft was inspired 
by similar efforts being done on a model rotor at 
the Rotorcraft Dynamics Laboratory 
(Longtrack) at Princeton's Forrestal campus l-J 
and aided by results reported from an 
instrumented AH-1G Cobra helicopter test by 
the NASA [3]. The design goals of the 
instrumentation system, and the impact ot 
flight safety and testing procedure on the 
realization of its mechanical and electncal 
components is described below. 


Rarir Instrumentation Goals 

Data Acquisition 

In order to provide a basis for comparison 
with model results, blade mounted 
accelerometers and strain gauges were s^ctM 
for the autogiro rotor experiments described. 
These sensors typically provide differentia 
outputs that require some form of amplification 
prior to sampling for data storage. On the 
model tests, both the sensor excitation voltages 
and differential outputs are transferred from 
the fixed to rotating frames using a slip nng 
assembly. Such a technique requires from two 
to four rings per sensor, plus the associated hub 
attachment for the sliprings and sensor wires. 
This type of installation was not possible, due 
to the inherent simplicity of the autogiro s 
’ hub. Since no torque is provided to the rotor, 
the primary aircraft control is performed 
through direct shaft tilt of the rotor, which is 
mounted in a pillow block attached to a sealed 
bearing assembly. Such an arrangement makes 
the main rotor shaft inaccessible for slipring 
attachment or routing of wires through its 
interior. Thus, the decision was made at an 
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early stage to telemeter the data from the 
rotating to fixed frames without the direct use 
of rotating brushes or slipring pick-ups. Use of 
a rotor-mounted telemetering system thus 
requires a co-located electronic power source, 
since no direct wire transfer was then possible 
from the aircraft fuselage. In order to reduce 
battery size on the hub, all integrated circuits 
used CMOS chips wherever possible, with 
power supplied from two standard 9-volt 
transistor batteries. A photograph of t 
instrumentation assembly is shown in Figure 2. 

An additional system requirement was 
that the instrumentation system not adverse y 
affect the mechanical integrity or 
aerodynamic performance of the rotor blades. 
The blades on the autogiro are stock Bensen 
factory-built blades constructed from aluirunum 
sheeting riveted to a solid spar. The airfoil is 
a Bensen design, having a flat underside and 
slight reflexed trailing edge, originally 
developed for construction by the homebuilder 
from plywood sheeting. Since access to the 
blade’s interior was not possible, multi- 
conductor ribbon cable for sensor signal routing 
was secured to the flat underside of the blade 
using a combination of doublesided sticky tape 
and epoxy, covered with one spanwise and 
several chordwise segments of mylar adhesive 
tape- 

Aerodynamic performance considerations 
dictated the scheme used to condition the 
sensor signals. Unlike the system used in 
Reference [3], cost considerations did not allow 
development of custom millimeter-thick 
integrated circuitry for sensor signal 
conditioning and amplification. Since stock 
integrated circuits (IC's) would be used, a 
minimum number of chips could be tolerated at 
each sensor station in order to reduce the 
adverse drag penalty from surface 
irregularities they introduce. Although sensor 
noise would be lowest for co-located sensors, 
amplifiers and analog to digital (A/D) 
converters, the associated multiplexing 
necessary would add at least an additional 
chip, bringing the total to three ICs at each 
sensor’s spanwise location. Thus, only the 
sensor's amplifier is located on the blade span, 
with the A/D and multiplexing operations 
performed by a single IC at the rotor hub. This 
arrangement also required only four wires to 
extend over the entire rotor radius. An 
installed accelerometer and amplifier are 
shown in Figure 3. 


Sensor signal A/D sampling rate was 
traded off against multiplexing capability for 
the digitizing of the rotor data, with the final 
choice using the TLC541 LinCMOS chip from 
Texas Instruments as the primary workhorse 
for conversion of the data to digital fonn. This 
IC provides 11 channels of multiplexing 
capability into an 8-bit A/D converter, with an 
equivalent throughput rate of 1,024 samples 
per second on each channel. Since the nominal 
rotor speed of the autogiro is 375 rpm, the 
system’s Nyquist rate (maximum digital 
bandwidth) is 82/rev, well beyond any 
potential dynamic or aeroelastic phenomenon 
one might expect to observe. For this reason, 
there are no anti-aliasing filters used prior to 
A/D conversion. 

The multiplexer chip, originally designed 
for interfacing with a microprocessor, allows 
for direct control of both channel addressing 
and serial bit output rate through reasonably 
complex interfacing with timing and address 
pins. In order to avoid the requirement of co- 
locating a dedicated microprocessor, the 
complex timing patterns and address lines were 
stored ("burned") into an EPROM, driven by a 
counter and dock. The individual data bits for 
each of the 11 channels are sequentially 
loaded onto the serial digital data bus, 
followed by a test channel for receiver 
synchronization purposes. This "synch" word 
is modified to produce a string of 9 bits that 
cannot be duplicated by any combination of the 
11 channel’s data words, thus providing a 
unique marker for each "frame" of data. 

To save on the total number of chips at the 
rotor hub, the counter that drives the EPROM 
interface chip does not reset to a predetermined 
value, and thus the first 1024 8-bit words are 
cycled over for each frame. Since the timing 
pulses only occupy the first two-thirds of this 
address space, a blank sector is available in 
the current frame of serial data for future 
expansion. This might consist of either 
additional rotating frame channels, or inter- 
woven fuselage sensor data. This latter 
method would require a serial decoder and 
synchronizer in the non-rotating frame on the 
aircraft, a method that was not used on the 
current design concept. A diagram of the main 
circuit functional blocks appears in Figure 4. 


Data Transmission a nd Storage 

In order to transmit the digital data from 
the rotor to the fuselage frame of the autogiro. 
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the serial bit stream was Manchester encoded 
so as to provide a timing reference for 
individual bit transitions. While this format 
is particularly suited for direct radio frequency 
modulation, the additional weight and power 
requirements of a radio transmitter were 
deemed unacceptable. Instead, the coded 
signal is fed directly into the rotor blade s 
aluminum structure, with the pickup signal in 
the fuselage frame merely consisting of a wire 
secured to the metal airframe. Thus, the entire 
aircraft is electronically isolated from the 
sensor power signals, with the rotor data 
transmitted directly through the main rotor 
shaft bearing assembly. Alternate schemes 
using infa-red diodes, while conceptually 
feasible, were not nearly as simple as this 
technique. 

Due to the limited amount of both ground 
support personnel and computational facilities, 
direct re-broadcast of the data from the 
autogiro to a ground station was not possible. 
Instead, appropriate signal modifications were 
made to the serial coded data so that it could 
be directly recorded onto an 8mm format video 
camera/recorder. The camera is a self- 
contained unit running from its own battery 
source, and is shock mounted to the autogiro 
directly in from of the pilot’s seat, to allow for 
convenient access to both the record buttons and 
tape eject mechanism. Despite the high 
frequency of the Manchester encoded serial bit 
stream (13 MHz), such signals are well within 
the bandwidth of conventional NTSC video 
standards. 

Post-flight data processing consist of 
playing back the recorded signals into a serial- 
to-parallel digital data synchronizer. The 
circuitry used to perform the data extraction, 
shown in functional form in Figure 5, consists of 
signal conditioning to standard logic levels, bit 
synchronization using a phase-locked loop 
(PLL), and data decoding with clocked flip- 
flop circuits. In order to provide for 
discrimination between channels, a 
search /synchronizer scheme using another 
EPROM and comparitor was employed. Since 
the frame word consisted of a unique 9-bit 
digital bit pattern, a shift register and bit 
comparitor were used to detect a ■’match" with 
this word indicating the beginning of the serial 
data frame. When the match was detected, 
the EPROM clock/counters were reset, and each 
channel was clocked through the shift register, 
combined with the four bits representing the 
channel count, and sent into a high-speed 


parallel digital data port on an IBM-PC/ AT. 
This data was then stored onto floppy 
diskettes for additional analysis and 
processing. After the EPROM cycled through 
its 11 channels of serial data timing, it entered 
a "re-match" state in which it looked for 
another frame start pattern. If this pattern 
was not found, the synchronizer would enter a 
"search" mode and warn the user (via light 
and buzzer) that synchronization was lost and 
data is invalid. If the pattern was matched, 
indications of a "locked" state would be given, 
and the process would cycle over each 
additional frame of serial data. 


/UiYiliarv Data Syst em Components 

In order to interpret the rotor structural 
dynamic data, some measurement of the 
operating state of the autogiro flight condition 
was required. Of particular importance is the 
rotor advance ratio, and the associated angle 
of attack of the rotor blades. As the detailed 
flow of the rotor is unavailable, several 
fuselage-mounted sensors were used to infer 
this information. These sensor data were 
routed to a pulse-width modulator (PWM) box 
having a capacity for 43 channels of analog 
data sampled at 20Hz per channel. This data 
acquisition system was previously used in 
flying qualities experiments at Princeton on 
variable-stability general aviation aircraft 
and gliders [4j. The PWM signal was in turn 
fed into the separate audio input channel of 
the 8mm video camera, thus providing 
separate but synchronous data recordings for 
both rotor and fuselage sensor data. 

The rotor speed was measured using a 
magnet mounted on the rotor hub, and a Hall- 
effect sensor on the rotor mast just below the 
hub pillow block assembly. The pulses were 
fed into a PLL circuit that output an analog 
voltage proportional to pulse frequency, with 
this information displayed to the pilot and 
sent to the fixed-frame data commutator 
system. Autogiro angle of attack and sideslip 
were measured through two vanes mounted on 
low friction potentiometers, and airspeed was 
taken from a cup anemometer. Rotor shaft tilt 
was measured from stick potentiometers, and 
rudder pedals were instrumented as well. A 
picture of the data system components is in 
Figure 6. 

Power for the PWM commutator was 
provided from four 9-volt batteries tied in 
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series, regulated to a nominal 28-volt DC level. 
Although larger power rechargeable batteries 
were available, their significant increase in 
weight deemed them unacceptable for initial 
flight tests and data system check-outs. 
Additional instrumentation, such as use o 
Princeton’s Inertial Measurement Unit 14], will 
require upgrading the available power source 
on the autogiro airframe. 

Since monitoring of the video tape index 
during flight could jeopardize safety through 
increased pilot workload, the FM 
communicators used between pilot and tow 
vehicle personnel was coupled into the audio 
input channel of both the on-board video 
recorder and the tow car’s video camera. Since 
this same space on the tape is used by the 
fixed-frame sensors, voice markers were stored 
only during push-to-talk operations of the 
pilot. Simultaneous voice recording allowed 
for approximate synchronization between the 
signals from both cameras during postflight 
analysis. 


preliminar y Testing 

As a means of both adjusting amplifier 
offsets and confirming data transmission, 
storage and retrieval/synchronization, an 

impact test was made of the rotor blade. An 
instrumented force hammer in conjunction with 
the blade-mounted accelerometers were 
recorded for a series of impacts at several 
spanwise locations in order to identify the non- 
rotating mode shapes of the autogiro rotor 
blade. Such information is essential for 
accurate post-processing of the accelerometer 
data, as is pointed out in 12]. 

Since the spanwise accelerometers are 
oriented to measure out-of-plane accelerations, 
as the blade deflects out of the plane of 
rotation, these sensors will measure both blade 
vertical acceleration as well as a component of 
rotational acceleration proportional to the 
local slope at the sensor's spanwise location 
(see Figure 7). This information can then be 
used in a processing scheme called a Kinematic 
Observer [6] to reconstruct rotor blade state 
variables. On the autogiro, however, the ratio 
between the measurement of fundamental 
teetering (flapping) acceleration and teetering 
displacement is constant for any spanwise 
station. For any accelerometer at spanwise 
location r, it will sense contributions from each 


mode’s acceleration <g(t)) and each mode’s 
displacement (g(t)) according to: 

«o 3rvj(r) 2 

accel(r,t) = gjft) + — r « gj(t) 

where ^(r) represents a particular blade 
natural mode shape. For the case of the rigid 
teetering mode, Tijto = r, giving contributions 

from the teetering mode Pj(t) as: 
rfr(t) + rfl 2 Pi(t) 

and for simple harmonic flapping at 1/rev, 
these two terms cancel. For this reason, a 
potentiometer was installed in the rotor hub to 
measure rigid teetering motion, and the 
accelerometers were positioned so as to 
maximize their sensing of the various higher 
blade vibratory modes. 

Flight operations for testing the autogiro 
instrumentation are currently underway at the 
inactive runway at Princeton University s 
Forrestal Campus. Since the autogiro has no 
engine, it is towed by a nylon rope attached to 
its nose from an automobile. While this results 
in limited continuous flight time, the runway's 
3000 foot length allows flight experiments of 
approximately 45 to 60 seconds duration, 
depending on wind direction along the tow 
path. Communications are kept with the tow 
vehicle driver and observer using FM 
transceivers, and a glider tow hitch with 
release may be operated by the pilot in the 
event of fouling of the tow line. 

Rotor data power has been tested to 
provide consistent data for over two hours of 
operation, and fuselage batteries have a 
roughly equivalent life. Battery life on the 
video camera used for serial data storage is 
slightly under an hour, resulting in 
approximately fifteen flights during a 
standard sequence of runs. Such a record 
. provides a wide range of test points for 
analysis, the results of which will appear in a 
forthcoming paper. 
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Figure 1: Instrumented Bensen autogiro (glider) 
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Figure 2: Digital data acquisition and telemetry circuitry 



Figure 3: Accelerometer and amplifier installation on blade underside 






(to video input) 


Figure 4: Rotor data acquisition circuit functional diagram 
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Figure 5: Digital data synchronization and parallel output functional diagram 














rc 6: Fuselage mounted airdata sensors and camera attachment 
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Fig. 7: Out-of-plane accelerometer sensing schematic 
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